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PART  ONE 

SOLUTION  PROCEDURE  AND  NUMERICAL  RESULTS 

I.  INTRODUCTION 

The  purpose  of  the  present  report  is  to  develop  an  H-field  solution 
which  can  be  used  alone  or  used  in  conjunction  with  the  E-field  solution 
in  [1]  to  obtain  a  combined  field  solution  which  is  hopefully  more  accurate 
than  the  one  in  [2].  Here,  H-field  solution  means  a  solution  to  the  mag¬ 
netic  field  integral  equation  for  a  perfectly  conducting  body  immersed  in 
an  incident  electric  field.  Similarly,  E-field  solution  means  a  solution 
to  the  electric  field  integral  equation.  In  [1]  and  in  the  present  report, 
solution  is  by  the  method  of  moments  with  a  Fourier  series  in  <j).  For  com¬ 
patibility,  the  expansion  functions  used  in  the  present  report  must  te  the 
same  as  those  in  [1],  In  [1]  and  in  the  present  report,  the  t  dependence 
of  the  expansion  functions  is  subsectional.  Pulses  are  used  for  the  (j)  com¬ 
ponent  of  the  unknown  electric  current  induced  on  the  surface  S  of  the  body 
of  revolution.  Triangles  divided  by  the  cylindrical  coordinate  radius  are 
used  for  the  t  component.  Here,  t  and  cj>  are  orthogonal  coordinates  on  S,  t 
being  the  arc  length  along  the  generating  curve  of  S  and  0  the  azimuthal 
angle. 

Computations  indicate  that  the  E-field  solution  of  [1]  is  more 
accurate  than  the  E-field  solution  of  [2]  and  that  the  present  H-field  solu¬ 
tion  is  nearly  as  accurate  as  the  H-field  solution  of  [2],  Hence,  the  com¬ 
bined  field  solution  obtained  by  putting  together  the  present  H-field  solu¬ 
tion  and  the  E-field  solution  of  [1]  should  compare  favorably  with  the 
combined  field  solution  of  [2]. 

Early  work  on  H-field  solution?,  summarized  in  [3] .  Two  recent 
H-field  solutions  are  presented  in  [4]  and  [2].  The  solution  in  [A]  uses 
pulse  expansion  functions  for  both  t  and  (}>  components  of  the  unknown  elec¬ 
tric  current  induced  on  S.  H-field  formulations  are  basic  ingredients  in 
solutions  to  problems  of  scattering  by  a  dielectric  body  of  revolution  [5], 
[6],  [7]  and  a  dielectrically  clad  conducting  body  of  revolution  [8].  The 
expansion  functions  in  [5]  and  [8]  are  triangles  divided  by  the  cylindrical 
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coordinate  radius  for  both  components  of  the  current.  The  H-field  formu¬ 
lation  of  [2]  is  used  in  [6]  and  [93 •  The  expansion  functions  in  [7]  are 
pulses  for  the  <j>  component  of  current  and  shifted  pulses  divided  by  the 
cylindrical  coordinate  radius  for  the  t  component. 

The  magnetic  field  integral  equation  for  a  perfectly  conducting 
body  of  revolution  illuminated  by  an  incident  magnetic  field  H_  is 

-n  x  hS(J)  -  n  x  r1  just  inside  S  (1) 

where  S  is  the  surface  of  the  body  of  revolution  and  n  is  the  unit  vector 

which  is  normal  to  S  and  which  points  outward  from  S.  In  (1),  £  is  the 

s 

electric  current  induced  on  S  and  H  (J)  is  the  magnetic  field  due  to  J 
radiating  in  the  absence  of  the  body  of  revolution.  It  is  assumed  that 
the  medium  outside  S  is  linear,  isotropic,  and  homogeneous.  "Outside  S" 
means  outside  the  body  of  revolution  and  "inside  S"  means  inside  it. 
"Absence  of  the  body  of  revolution"  means  that  the  body  of  revolution  has 
been  removed  and  that  the  ensuing  void  has  been  filled  with  the  outside 
medium.  The  incident  magnetic  field  H*  is  the  field  which  would  exist  in 
the  absence  of  the  body  of  revolution. 

Equation  (1)  states  that  the  tangential  components  of  the  total 
s  i 

magnetic  field  (H  +  H  )  vanish  just  Inside  S.  Consequently,  (1)  is 
valid  only  if  S  is  closed.  By  contrast,  the  E-field  integral  equation 
(40)  of  (2]  is  valid  for  both  open  and  closed  surfaces  of  revolution. 


II,  METHOD  OF  MOMENTS  SOLUTION 

Equation  (1)  is  solved  by  the  method  of  moments.  The  moments 
solution  to  (1)  is  approached  by  writing 


J  =  l  (  l  It,  J* 

-  £  J  nj  -nj 


J“.  +  l  I 


<t>  j<t>  - 

nj  -nj 


(2) 


where  and  are  expansion  functions  defined  by 
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u. 


Ti(t) 

P 


ejn<j) 


P.(t) 

_ M. _ - 


ej# 


fn  =  0 ,  +1 ,  + 2 , . . . 

1 1  -  1.2,...  P-2 

(n  =  0 ,  +1 ,  +2 , . . . 
j  =  1,  2,  ...  P-1 


(3a) 


(3b) 


and  where  I*",  and  1^,  are  unknown  constants  to  be  determined.  The 
nj  nj 

subscript  j  which  runs  from  1  to  either  (P-2)  or  (P-1)  in  (3)  is  not  to 
be  confused  with  the  j  which  appears  in  the  argument  of  the  exponential 
in  (3) .  The  latter  j  is  /-T.  The  variables  t  and  <J)  are  orthogonal 
curvilinear  coordinates  on  S.  Specifically,  t  is  the  arc  length  along 
the  generating  curve  of  S  and  cf>  is  the  azimuthal  angle.  The  generating 
curve  of  S  is  the  plane  curve  which,  when  rotated  about  the  z  axis, 
generates  S.  In  (3),  u^  and  are  unit  vectors  in  the  t  and  <{)  directions, 
respectively.  These  directions  are  chosen  so  that 


n  =  u, 

-  -Cj) 


X  ^t 


(4) 


The  function  T^  (t)  is  the  triangle  function  defined  by 


f~(t  -  t~)/A 


j 


T;(t)  »  <  (t4J_0-t)/A 


j+2  "  j+1  * 

0 


1  C  1  Cj+1 


Vi 1 c 1  V2 


otherwise 


(5) 


and  P^(t)  is  the  pulse  function  defined  by 


PjCO 


1  *  h-^Vi 


(6) 


0  , 


otherwise 


Here,  t~  ,  t2>  tp  are  consecutive  but  not  necessarily  equally  spaced 
points  on  the  generating  curve  of  S.  t^  is  the  value  of  t  at  the  be¬ 
ginning  of  the  generating  curve  and  tp  is  the  value  of  t  at  the  end. 
Henceforth,  the  generating  curve  is  assumed  to  be  a  series  of  straight 
line  segments  joining  the  points  t^,  t^,...  tp.  In  (5), 
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Aj '  Vi '  b 


(7) 


In  (3),  p  is  the  distance  of  the  point  t  on  the  generating  curve  from 
the  z  axis  and  is  the  value  of  p  at  t  =  t^  where 

Cj  '  2  (tj  +  VP  (8) 

It  is  assumed  that  all  the  Pj’s  are  positive. 

The  moments  solution  to  (1)  is  obtained  by  substituting  (2)  into 
(1) ,  integrating  over  S  the  dot  product  of  (1)  with  each  member  of  a 


set  of  testing  functions,  and  solving  the  resulting  matrix  equation  for 

Ifc.  and  I* 

*3  nj 

defined  by 


It.  and  1^,.  There  are  two  kinds  of  testing  functions,  W*".  and  . 
nj  nj  &  -ni  ~ni 


t  Ti(t)  -1n^ 


J  pi(t)  -an<j> 

VT  4  =  u,  — -  e  J 


•nai  ^  P4 


(9a) 


(9b) 


The  matrix  equation  for  1^  and  1^  decomposes  into 


rtt 

Yt<r 

r?r 

-it 

n 

n 

n 

n 

<j>t 

t* 

|i<j> 

n 

n  _ 

L*- 

_  n  _ 

n  =  0,  +1,  +2, . . . 


(10) 


where  the  superscripted  I  's  are  column  vectors  and  the  superscripted  Y  's 

n  ->-t  t  n 

are  submatrices.  The  jth  element  of  1^  is  1^  and  the  jth  element  of 

1^  is  T^,.  The  ith  elements  of  and  I ^  are  given  by 
n  nj  n  n 


T 

1  •  •  .•TT> 


2tt 


l1; 

ai 


dt  T±(t) 


d$(u^  *  Hi)e_:5n<^ 


(Ha) 


2-rr 


r*  -  -  -L 


^i-"^  I  dtpPi(£) 


d<Ku.  •  H*)e 
t  ~~ 


(Hb) 


The  second  subscript  on  I  in  (11)  denotes  the  ith  element.  The 
superscript  i  on  I  in  (11)  is  the  same  as  that  in  (10).  This  super¬ 
script  indicates  dependence  on  the  incident  magnetic  field  Hi. 

For  a  particular  value  of  n  in  (10) ,  the  matrix  of  the  super¬ 
scripted  Y  's  on  the  left-hand  side  of  (10)  is  called  the  moment 
n 

matrix  for  that  value  of  n.  The  ijth  elements  of  the  superscripted 
Y^'s  are  given  by 


'i+2 


=  * 


T±(t)  Tj (t) 


i+2 


dt  +  k 


1+2 


i  _ 


dt  T^t) 


dt’  T  (t') 


[((p,-p)cos  v'  -  (z'-z)sin  v')  0  cos  v’j 


(12a) 


~  i+1 

„  IT 


t . 


(yy  ) 

Un  Jij  p 


dt  p  P±(t) 


+2 

dt’  T.(t')(p'sin  v  cos  v’ 
i-  1 


'i  j 

-  p  sin  v’  cos  v  -  (z’-z)  sin  v  sin  v’)G 


_  £i+2 

» f-  f 

n  ij  Pj 


j+1 


dt  TA(t) 


dt’  p’  P^ (t ’ ) (z’-z)G3 


(12b) 


(12c) 


(Y**) 


i+1 


n  ij 


'i+1 


PPl(t)Pj(t)dt+—  J 


dtpP^t) 


+1 


dt’p’P  (t1) 


[((p’-p)cos  v  -  (z’-z)sin  v)G^  +  G^p’  cos  v] 


(12d) 
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The  derivation  of  (12)  is  not  included  here  because  it  is 
that  of  (20)— (23)  of  [2].  Actually,  (12)  was  obtained  by 
f^(t)  by  —  T^(t)  in  (20)  and  (22)  and  by  P^t)  in  (21) 

and  by  replacing  f.(t')  by  -4-  T.(t')  in  (20)  and  (21)  and 
J  ^  J 

in  (22)  and  (23).  The  equations  (20)-(23)  referred  to  in 
sentence  are  in  [2]. 


similar  to 
replacing 
and  (23) 


by^h(t’> 

the  previous 


In  (12),  k  is  the  propagation  constant  in  the  medium  outside  S. 
Also,  v  is  the  angle  between  the  tangent  to  the  generating  curve  and 
the  z  axis,  v  is  positive  if  the  generating  curve  is  moving  away  from 
the  z  axis.  Otherwise,  v  £  0.  In  (12),  p,  z,  and  v  are  evaluated  at 
the  point  t  on  the  generating  curve  and  p',  z’,  and  v'  are  the  values 
of  p,  z,  and  v  at  the  point  t'.  The  subscripted  G’s  in  (12)  are  given 
by 


Gl  =  2 


f 


G  sin  (<f>/2)  cos(n<j>)d<j> 


(13a) 


G2  = 


G  cos  <1>  cos(n<|>)d<j> 


(13b) 


G3 


G  sin  <J>  sin(n<j>)d<j> 


(13c) 


where 


where 


G  = 


-  1  +  .IkR  “jkR 

3  3  e 

k  K 


R  = 


(p'-p)2  +  (z'-z)2  +  4pp'sin2(())/2) 


(14) 


(15) 


The  regions  of  integration  in  (12a) -(12c)  overlap.  To  avoid 

integrating  more  than  once  over  the  same  region,  we  account  for  the 

contributions  in  (12a) -(12c)  by  regions  of  integration  rather  than  by 

matrix  elements.  If  A  is  the  region  which  t  <  t  <  t  ,, ,  then  A 

q  6  q-  -  q+1  q 

contributes 
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_  f  T/t> 

(Yct)  =  n  _± - J - dt 

n  'ij  p 


to  the  single  integral  with  respect  to  t  in  (12a) .  Because  the 
domains  of  the  triangle  functions  are  limited,  (16)  is  nonzero 
only  for 


i  =  q-1,  q 


i  j  0 
i  i  P-1 

j  q-i.  q 
j  ^  0 
j  7*  P  -1 


The  dot  on  the  left-hand  side  of  (16)  indicates  that  (16)  is  not  all 
of  (12a)  but  only  the  contribution  due  to  A  . 

q 

If  A  is  the  region  for  which 

pq 

t‘  <  t  c  t" 

p  -  ~  p+1 

t“  <  t*  <  t“  , 
q  -  ~  q+1 

then  the  contributions  to  (12a)-(12c)  due  to  A  are  given  by 


(17a) 


(17b) 


(Y“)±j  -  k3  j  dt  T^t)  j  dt’T^t')  [((p'-p)cos  v'  -  (z’-z)sin  v’)G2 


-  p  cos  v'} 


(Y^)  ,  =  dtpP  (t)  dt’  T.(t,)(p'  sin  v  cos  v' 

n  pJ  P„  J _  p  J_  i 

F  t  t 

p  q 


(18a) 


-  p  sin  v'  cos  v  -  (z'-z)sin  v  sin  v')G 


(18b) 


8 


P+1 


q+1 


^t4>>4n=4r-  dtT,(t)  dt’p»  p  (t’)(z'-z)G. 
n  iq  p„  i  q  i 

s  t  t 

p  q 


where 


and 


(18c) 


(19a) 


(19b) 


The  asterisks  on  the  left-hand  sides  of  (18a)-(18c)  indicate  contri¬ 
butions  due  to  Apq.  For  instance,  (18a)  is  not  all  of  (12a)  but  only 
the  contribution  due  to  A  .  Equation  (12d)  is  suitable  for  calcula- 

pq 

tion  as  is  because  it  has  no  overlapping  intervals.  However,  to  ob¬ 
tain  a  notation  consistent  with  that  in  (18a) -(18c),  we  replace  ij  by 
pq  in  (12d) .  The  result  is 


(Y#) 

n  pq 


’p+1 


2  k 

pP  (t)dt  +  — 
P  P„P 


p  q 


j  dtpPp(t) 

t~ 

p 


J  dt,p,pq(t') 

t” 

q 


where 


[((p'-p)cos  v  -  (z’-z)sin  v)G£  +  G^p!  cos  v] 


The  fact  that  the  pulses  (6)  do  not  overlap  was  used  to  obtain  the 
S  term  in  (18d) . 

pq 


(18d) 


(20) 


V.t '  V  -*•  ,  M' '  '* »  ^  '  '  ‘  1  V  '  " ' 


Each  integral  with  respect  to  t  in  (18)  is  approximated  by 

sampling  the  pertinent  integrand  at  t  =  t  and  multiplying  by  Ap. 

Next,  the  definitions  (5)  and  (6)  of  the  triangle  and  pulse  functions 

are  substituted  into  (18).  Because  the  portion  of  the  generating 

curve  for  t  <  t*  <  t  ,,  is  assumed  to  be  a  straight  line,  v*  is  con- 
q  -  -  q+1 

stant  there  so  that 


p*  =  p  +  (t*  -  t  )sin  v 
q  q  q 


t  <  t'  <  t  .. 
q  -  -  q+1 


(21) 


z'  =  z  +  (t'  -  t  )cos  V 

q  q  q 


where  z  is  the  value  of  z'  at  t*  =  t  .  Equations  (21)  are  also  sub- 

q  q 

stituted  into  (18) .  As  a  result,  (18)  becomes 


($“) 


k3A  Vl 


n  'ij  4  J 


„  ,  2(t’-t  ) 

»q-j _ . _ a 


(1  +  (-l)q“J  - — -q-)[((pq~Pp)cos  vq  -  (Zq-Zp)sin  vq)G2 


-  G-p  cos  v  ]  d  t  ’ 

1  P  q 


(22a) 


1k3A  q+1 


n  Pi  2  J 


(1  +  (-l)q-^  _ta^ 

'  '  - 7 — )  (p  sin  v  cos  v  -  p  sin  v  cos  v 

A  '  N  l  -rv  n  n 


q  p  q  p  q  p 

(22b) 


(z  -z  )sin  v  sin  v  )G.dt' 
q  P  P  q  3 


,  jk3A  (t'-t  )sin  v 

<*■  +  - h - +  (t'"tq)cos  vq)G3dt’  (22c) 


n  '  iq  2 

t 

q 


q  p 


trA 


'q+1 


(YW)  =  — 3-  6  +  k  A 

x  '  r\i*t  A  nn  r 


n  pq  pa  pq  p  J_ 


(t’-t  )sin  v 

(1  + - a- - -%(Pq-pp)cos  vp-;Zq-Zp)sin  vp 


(t'-t  )sin  v 

+  (t'-t  )  (sin  v  cos  v  -  cos  v  sin  v  ))G.  +  GnP  (1+ - - 3-)cos  v  ldt 
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The  ranges  of  values  of  i  and  j  in  (22)  are  given  by  (19) .  Because  p 

is  linear  in  t,  the  process  of  sampling  the  integrand  at  t  =  t  and 

P 

multiplying  by  Ap  gave  the  exact  value  of  the  first  integral  in  (18d) . 

The  integral  in  (16)  can  be  approximated  either  by  sampling 

the  integrand  at  t  =  t  and  multiplying  by  or  by  applying  Gaussian 

quadrature.  Gaussian  quadrature  was  chosen  because  it  gave  more 

accurate  values  of  electric  currents  induced  on  the  sphere  and  cone- 

sphere  examples  in  Section  V.  Thanks  to  (5)  and  (21),  the  Gaussian 

quadrature  formula  for  the  integral  in  (16)  becomes 

„  (n  )  00  9 

TfA  t  A?  C  (1  -  x,  ) 

fYtt)  =— SL  V  J: - 2 - 

Un  ;q-l,q-l  8  ^  0£ 

_  (n. )  (n  )  (n  ) 

ttA  t  A,  (1  -  x0  )  (1  +  x„  ) 

/*ttv  =  fYtS  =  — 1  y  - - K - - — 

Un  q,q-l  Un  ;q-l,q  8  ^ 


(23a) 


(23b) 


(n  )  (n  )  , 

ttA  t  A.  t  (1  +  x?  ) 

(Yct)  y  -i - - 

Un  q,  q  8  ^ 


(23c) 


where  |30  =  P  + 

A/  (J 


~2  Sln  Vq 


(23d) 


(n  )  (n  ) 

The  abscissas  x£  and  weights  A£  in  (23)  are  given  in  Appendix  A 
of  [10]  for  several  values  of  nfc. 

Equations  (22)  can  be  rewritten  as 

t((Pq-pp)c°s  vq  -  (yzp)sin  »q)021  -  Sy  Pp  ™  »,! 


+  (-1) 


(Y**)  „ 

n  'pj 


k3A  A 


jk3A  A 


[((Pq~Pp)C0S  Vq  "  (YZp)sin  Vq)G22"G12PpCOS  Vq]  ^ 


, k  3  (n  sin  v  cos  v  -  p  sin  v  cos  v  - 
4  q  p  q  p  q  p 


-  (z  -z  )sin  v  sin  v  )(G„-  +  (-1)^  3G„0) 
q  p  p  q  31  3z 


(24b) 
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ik3A  A  A  sin  v  A  cos  v  A  sin  Vo 

t  <  VV  <=3!  +  G3^  +  (G32  +  -^Tof-  V  > 


(Y^). 
n  iq 


(24c) 


nA  b.  u  u  ~0  n 

cd*)„  -  x*  «  „  +  — K(p -pJcos  V(V.l8f“  V<021+  2o  G22) 


k3A  A 


A  sin  v 


‘n  '  pq  Pp  pq  2 


q  p  P  q  P 


A  sir.  v 

+  (sin  v  cos  v  -  cos  v  sin  v_)(Go0  +  q  G^)  + 


A 

_c 

2 


q  p 


A  sin  v 


,  oxn  *  yv»22  •  2p  23 

H  q 

A  sin  v 


+  pqcos  vp(Gu  +  "q  1  Gl2  +  (  £  2p“£)  G13)] 


(24d) 


where 


3+1 


r  =  Aa 

ma  A 

q 


(t’-t  )a_1G  (t’-t )dt’ 

q  m  q 


m  =  1,2,3 
a  «  1,2,3 


(25) 


Here, 


r 


Gx(u)  =  2 


5(u,<t>)sin  (<t>/2)  cos  (ncj>)d<j> 


(26a) 


°2<u)  ' 


G(u,<i>)cos  4>  cos  (nc|)) d4> 


(26b) 


■f 


G^(u)  -  j  G(u,<J))sin  <p  sin  (ruj>) dcf) 


where 


where 


G(u,«  -  e'3kS<UlW 

txM) 


(26c) 


(27) 


R(u,<j>) 


/ 


where 


(p’-pp)2  +  (z'-Zp)2  +  4ppp’  sin2(<t>/2) 


P*  =  p  +  u  sin  v 

h  q 

z 1  =  z  +  u  cos  V 


(28) 

(29a) 

(29b) 
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III.  NUMERICAL  EVALUATION  OF  THE  INTEGRALS  G 

- - - met 

The  integrals  Gma  of  (25)  are  now  evaluated  by  means  of  either  pure 
Gaussian  quadrature  or  Gaussian  quadrature  fortified  by  the  method  of 
elimination  of  the  singularity  of  the  integrand.  The  proximity  of  the 
singularity  of  G(u,<)>)  of  (26)  to  the  region  of  integration  falls  into 
one  of  four  different  cases  defined  in  this  section. 


If 


P  +  q 


-  c  A  <  d 
2  t  q  —  o 


c .  p 

<J>  q 


Case  1 

Pure  quadrature 


(30) 


then  pure  Gaussian  quadrature  is  used  to  evaluate  the  integrals  in  (25) 
and  (26).  In  (30),  cfc  and  c^  are  constants  for  which  the  \alues 


ct  =  2- 


c,  -  0.1 
<P 


are  suggested,  ^lso  in  (30), 

* 


d  -■ 
o 


i*i  1 

t  <  4  A 
1  o'  —  2  q 


(31) 


(32) 


yRFT 


2  *  2 

Aq) 1  +  (d  r 


i  *  i  i  . 

t  >  -  A 
o'  2  q 


where 


t  =  (p  -p  )sin  v  +  (z  -z  )cos  v 
o  q  p  q  q  p  q 


d  =  I (p  -p  )cos  v  -  (z  -z  )sin  v 
4  P  q  q  p  q 1 


(33) 

(34) 


The  quantity  d  defined  by  (32)  is  the  minimum  value  of  R(u,<£)  for  the 

ranges  |u|  <  -r  A  and  0  _<  <|)  <_  u  of  values  of  u  and  <j>  in  (26).  Equations 
l  q 

(30)-(34)  were  taken  from  pages  24  and  25  of  [1], 


The  Gaussian  quadrature  formulas  for  the  integrals  in  (25)  and 
(26)  are 

°t  (n  )  ’r.'  o-l  ,  (n  )  '  ’ "  *  1,2,3 

Gmct  "  ,?  ,  V  >  °m(2  Aq  V  >  *' 

11  '  a  -  1,2,3 

V. 

Il4>  (n.)  2  1 

G-^u)  =  it  l  A^  ™  G(ut<J>£>sin  (y  4y)cos  (n^) 

i  £=i  ~  ^  * 


(35) 


(36a) 


7T  *  (nd,) 

G0(u)  =  l  A0  ™  G(u,<j>Jcos  <k  cos  (aj)-) 


?  L  5, 

Z  £=1  * 


<J)  (il.  ) 

g3(u)  =  y  l  A^  ^  G(u,4>2)sin  4>&  sin  (n<{y) 
&=1 


where  G(u,<J)^)  is  given  by  (27)-(29)  with  4>  replaced  by  <jy.  Here, 

$0 


1  (x(V  +  d 


(nj 


(36b) 


(36c) 


(37) 


<V  <»t>  <V 

The  abscissas  x^,  and  x^  and  weights  A^t  and  A^  in  (35)  and 
(36)  are  given  in  Appendix  A  of  [10]  for  several  values  of  n^  and  n^. 

If  (30)  is  not  true,  then  the  singularity  due  to  the  kernel  (27) 
is  eliminated  [11]  before  applying  Gaussian  quadrature  to  the  integrals 
in  (25)  and  (26).  The  singularity  analysis  presented  here  is  more  gen¬ 
eral  than  that  in  [7]  because  the  one  in  [7]  was  carried  out  only  for 
p  =  q.  Substitution  of 

2  2 

-jkR  .  kR 

e  J  =  1  -  jkR - j — 

Into  (27)  reveals  that,  near  R  =  0,  (27)  behaves  like  a  function  GS 
given  by 

„s  _  1  +  1 


,  3n3  2kR 
k  R 


(38) 


Three  different  methods  of  eliminating  the  singularity  are  used. 
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p  i*  q 


—  c  A  >  d 
2  t  q  o 

c.p  <  d 

<t>  q  -  o 


Case  2 
Method  1 


(39) 


then  method  1  is  used.  In  method  1,  the  orders  of  integration  in  (25) 
are  interchanged  to  obtain 


G  =2 
la 


IT 

0 

,TT 


H^(<j>)sin  (tj>/2)cos  (ncf>)d<j> 


J2a 


G3a  " 


H^((j))cos  cos(n({>)d<f> 


0 


>  a  =  1,2,: 


fTT 


H  (cj))sin  4>  sin(n4>)d<j; 


a 


(40a) 


(40b) 


(40c) 


where 


2  a 

H  (<f>)  =  (~~) 
a  Y  A 

q 


a+1 


(t'-tq)a_1  G(t'-tq,<j>)dt’ 


(41) 


In  method  1,  the  singular  part  of  the  integrand  in  (41)  is  integrated 
analytically,  the  remaining  part  is  integrated  by  means  of  Gaussian 
quadrature,  and  then  Guassian  quadrature  is  applied  to  the  integrations 
with  respect  to  cf)  in  (40) .  '  The  singular  part  of  a  singular  integrand 
is  a  function  which  behaves  like  the  integrand  near  its  singularity. 
Evidently,  the  singular  part  is  not  unique. 


Guided  by  (38)  in  our  choice  of  the  singular  part  of  the  inte¬ 
grand  in  (41),  we  write 


q+1 


2  a 

B  ■.-)  -  (f) 

q 


(f-o 


)<x“1  (G(t'-t  ,(j>)  -  - -rr^T-r r ~  ■  > 

q  q  k3R3(t'-t  ,<J>)  2kR(t 


)dt'+i 


a 


(42) 
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where 


ln  =  (Hr. )a  I  (t'-t  )^1  (---.-■— - +?ki77»'-t  TPdt'  (43) 

a  Aq  J-  <  k3R3(t*-t  ,♦)  2kR(t 


Application  of  Gaussian  quadrature  to  the  first  integral  in  (42)  gives 


H  (<j>) 
a  T 


nt  (n  )  (n  )a-l  -  (n  )  , 

l  CV')  (G(|  A*/  ,«  -  — - (T)— 

k  R  ,4>) 


- — _ — .)  +  I 

1  (nt)  a 

2kR(|  A  x  C  ,<}>) 


Equation  (43)  can  be  rewritten  as 


vK 

Xa  =  (f)Ct  t~3,  2  2.3/2  +  ...  2  2.1/2,d" 

q  ■'1.  k  (w  +d  )  2k  (w  +d  ) 

t  -  ^rA 
o  2  q 


where 


t  =  (p  -p  )sin  v  +  (z  -z  )cos  v  +  2p  sin  v  sin  (<j>/2) 

o  v  q  p  q  q  p  q  p  q 


(46a) 


A1  1  , 2 
d  =  r  -  t 
pq  o 


(46b) 


r2  -  (P  -p  )2  +  (z  -z  )2  +  4PP  sin2(<}>/2) 
pq  q  p  q  P  P  q 

Application  of  formulas  200.01.,  200.03.,  and  sequel  in  [12]  to  (45) 


(46c) 


gives 


2  w  ,  1 


h  -  k r  7TT  +  2  l08("+r>  . 

q  k  d  r  _  1 


t  +  -k 

o  2  q 


(47a) 


t  -  xA 
o  2  q 
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2  2 


te  JL 

2  "  kr 


V  2^  q  2t 

- -  i 

A 

•  1  q 
t  -  ih 
o  2  q 


I,  = 


2  3 

q 


2  2 

^  +  (1  "  ^-)log(wfr) 


t  -HrA 
o  2  q 


4t  I0  2t  „ 

-  2-1  _  (^)2  ! 

A  v  A  ^  1 


o  2  q 


(47b) 


(47c) 


where 


=  yjv' 


2  2 
+  d 


(47d) 


If  4>  =  0  in  (42),  then  (30)  guarantees  that  the  variable  t* 
of  integration  in  (42)  passes  close  to  the  singularity  of  G(t'-t^,<j>) . 

In  this  case,  (44)  and  (47)  are  necessary.  However,  as  cj>  moves 
toward  ir,  r  ^  may  become  so  large  that  the  variable  t'  in  (42)  never 
comes  close  to  the  singularity  of  G(t'-t^,<j>) .  In  this  case,  (44) 
and  (47)  are  no  longer  necessary.  Their  calculation  is  always  time 
consuming  regardless  of  the  value  of  r  .  Even  worse,  the  calculation 
of  (44)  and  (47)  is  subject  to  excessive  roundoff  error  whenever  r  is 

pq 

appreciably  larger  than  A  .  It  was  decided  to  restrict  use  of  (44) 

q 

and  (47)  to 


r  <Jci  +h 
Pq  2  t  q  2  q 


(48) 


For  larger  values  of  r  ,  recourse  is  to  the  Gaussian  quadrature 

pq 

formula 


.  nt  (n  )  (n  )  a-1  .  (n  ) 

-  „L,  v  <*r  >  Vr 


(49) 


The  term  -^A  on  the  right-hand  side  of  (48)  is  necessary  to  assure 

that  the  distance  R(t'-t  ,<j>)  used  in  (41)  is  at  least  as  large  as 
1  q 

2  ctA  for  all  values  of  t*  in  (41)  before  (49)  is  invoked.  This  dis¬ 
tance  enters  (41)  through  (27). 
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0  <  b  :$  <T/l2  (54) 

The  left-hand  side  of  (54)  is  a  consequence  of  the  fact  that  |x|  >_  | sin  x| 

for  any  value  of  x.  The  right-hand  side  of  (54)  was  obtained  by  setting 
2  2 

a  =  PpP1'}*  in  (53b).  Substitution  of  (52)  into  (38)  and  subsequent  expan¬ 
sion  of  (38)  in  powers  of  b  reveals  that,  near  <}>  =  0,  GS  behaves  like  the 
ss 

function  G  given  by 

ss  1  -1  pnp,<J)4 

GSS  =  — - —  +  — =—  +  —2- -  C55'v 

.3  3  2ka  01  3  5 
k  a  8k  a 

Guided  by  (55)  in  our  choice  of  the  singular  parts  of  the  inte¬ 
grands  in  (26) ,  we  approach  method  2  by  rewriting  (26)  as 


G1(u)  =  2  (G(u,<j>)sin  (<j>/2)cos(n<|>) - y~_)dc{>  +  — d<j> 


(56a) 


7T 

G2(u)  -  j  (G(U,«C08  *  cos(n«  -fj  -  ^  "  ^Tr)d* 

JQ  k  a  2k  a  8k  a 


f?r  i  i  /  2.1n.2  p  p’<J> 

+  (  1  +  -A-  (n  +!)({>  PK  T  y, , 

\3  3  +  2ka  ,.3  3  a,3  5  )d<i> 

k  a  2k  a  8k  a 


(56b) 


fir  2  j-tt  2 

G„(u)  =  (G(u,<)>)sin  cj>  sin(n<J>)  -  )d(|>  +  dcj) 


(56c) 


There  are  two  integrals  on  the  right-hand  side  of  each  of  equations 
(56a),  (56b),  and  (56c).  If  Gaussian  quadrature  is  applied  to  the  first 
integral  and  if  the  second  integral  is  evaluated  analytically  by  using 
formulas  200.01.,  200.03.,  and  sequel  in  [12],  then  (56)  becomes 


G1(u) 


<f>  (n.) 

V  A  (P 


I  9  (G(u,<()£)sin2(-|  41^)008 (n^) - j-y) 

=1  Air  a 


3  " — 372  tr  “  lo8(w+r) ] 

2kJ(ppp')J/2  r 


(57a) 


'T" 
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7T  ^  (V  1  1  (n2+l)4>n  P  P'4 

g2(u)  =  2  i  h  (G(uA)C08  Vos(nV  "TnkT +  -nr—  “  frr) 

l-l  k  a^  £  2k  8k  a^ 

+  2  3  ^  3/2  +  +  3(Al>3/2  g  ~  ' 

A3(ppp’)3/  r  2k(ppp’)1/2  2k3(ppp')3/2 


TTTJn  [g(!+3”  ? 


(57b) 


Sk^CPpP*)-'-  3r 


n 


$  (y 


G3(u)  =|  £  (G(u,4>£)sin  4)^  Sin(n4>£) - y-j) - 3 - - — 3/2  ^  ”los(wfr)] 

£=1  k  a^  k  (p  p') 

(57c) 


where  a^  is  the  right-hand  side  of  (53a)  evaluated  at  4>  “  4>n-  In  (57), 


w 


Wp  p’ 
_  P 


/(pf-p  )2  +  (z'-z  V 


(58a) 


r  =  yl  +  w2 


(58b) 


/2  2 
(p'-p  )  +(z'-z  ) 

P  P 

of  R(u,4>)  with  respect  to  4>  in  (26)  satisfies 


(p'-pp)2  +  (z’-zp)2 


Vq 


(59) 


The  distance  R(u,4>)  enters  (26)  through  (27).  If  (59)  is  not  true,  then 
G  (u)  is  approximated  by  the  pure  Gaussian  quadrature  formulas  (36) . 

In  method  2,  G  is  obtained  by  substituting  either  (57)  or  (36)  into  (35). 


If 


p  =  q 


Case  4 
Method  3 


(60) 


then  method  3  is  used  to  obtain  G^.  Since  p=q,  only  G^,  G.^»  G^ 
and  G^^  are  needed  in  order  to  calculate  (24) .  Prom  (40)  and  (41) ,  these 
subscripted  G's  are  given  by 
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2  a  r  f^q  a-1  2 

^1  =  2^&~ )  d<|>  du  u  G(u,<jb)sin  (4>/2)cos  (txcj>)  s  a  =  1,2,3 

"  °  -k 

2  q 


(61a) 


1. 

2  a  r  f2q  a-1 

^3a  =  ^~E~^  d<j>  du  u  G(u,<|>)sin  (j>  sin(n<j>),  a  =  2,3 

q  L  1. 

0  “  7An 

2  q 


(61b) 


where  G(u,v)  is  given  by  (27)  with' 


2  2 
R(u,4>)  =  \J u  +  4pq(pq+u  sin  v^)sin  (4>/2) 


Substitution  of  (38)  for  G(u,(J>)  in  (61)  reveals  that,  in  (61),  the  only 
integrand  which  is  not  bounded  in  the  vicinity  of  u  =  4>  =  0  is  the  inte¬ 
grand  for  a  =  1  in  (61a).  In  method  3,  the  singular  part  of  this  inte¬ 
grand  is  integrated  analytically  with  respect  to  both  u  and  the  re¬ 
maining  part  is  integrated  by  means  of  method  1,  and  then  method  1  is 
applied  directly  to  the  integrals  (40)  and  (41)  for  G12>  G^j  G^,  and 


Equation  (61a)  for  G^  is  rewritten  as 


Gin  =  1  +  I, 
11  a  b 


where 


tv  wA 

1 2  q 


h-rj  , 

0  -  k 

2  q 


du.(G(u,<j>)sin  (4>/2)cos (n<j>)  - 


3.  2,2.2.3/2)  (64) 

4k  (u  +pqtj)  ) 


1. 

1  ^  f2  q  A2 

Tb  =  T  g  ^  du  ,3'  2  ,  2  2  3/2 

q  0  -  k  <U  +  Pq^  > 


It  is  evident  from  (38)  and  (62)  that  the  integrand  in  (64)  is  bounded. 
Application  of  formulas  200.03.  and  200.01.  of  [12]  to  (65)  gives 
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1  *  -o~T  log 

b  k3p3 


2p  IT  /  2Pfl7r  2 

+  V 1  +  (T1-)2 


Integration  with  respect  to  u  of  the  second  term  in  (64)  by 
means  of  formula  200.03.  of  [12]  gives 


(66) 


I  ='2 
a 


H1  (4>)sin  (4>/2)cos(n<J>)d$ - r-v 

1  k  D 

K  Pq  0 


d<{> 


(67) 


where  H^)  is  given  by  (41).  The  first  integral  in  (67)  is  the  right- 
hand  side  of  (40a)  for  a  =  1.  Application  of  method  1  to  this  integral 
and  to  those  in  (40)  and  (41)  for  G^,  Gl3>  G32,  and  G^,  evaluation  of 
the  second  integral  with  respect  to  <f>  in  (67)  by  means  of  Gaussian  quad¬ 
rature,  and  use  of  (63)  and  (66)  give 

(n.) 

n4>  (nA)  2  i  *  >  A/ 

G.  =  tt  l  A0  9  H.  (()>.)  sin  (^-(McosfaK) - srv  l  ,  — 

11  1  1  £  2  *  **  2kV  l- 1  rZ  A~7 

q  V  »2+e> 


.3  3 

kPq 


log 


2p  IT  /  2p  IT  ' 


(68a) 


U  »  ^  ^ 

Id  ■”  /,  Al  41  W8ln‘(2  't,i)c°8(“V 

jo25! 


a  -  2,3 


<o 


(68b) 


(68c) 


G3a  “  f  ^  H  *  Wsin  *2 

where  H  (<j>D)  is  given  by  either  (44)  or  (49)  depending  on  whether 

Ot  X* 

(48)  is  true. 

Equations  (68b)  and  (68c)  and  the  first  sum  in  (68a)  were  taken 
directly  from  (50).  The  rest  of  the  right-hand  side  of  (68a)  is  due  to 
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the  manipulations  (63) — (67) .  The  numerical  integration  with  respect 
to  <j>  required  to  obtain  (68a)  is  of  the  bounded  integrand  in  (64) 
rather  than  of  the  unbounded  integrand  in  (40a).  Hence,  the  portion 
of  (68a)  in  excess  of  (50a)  fortifies  the  numerical  integration  with 
respect  to  <J>.  Having  been  obtained  through  analytical  integration 
with  respect  to  u,  this  portion  of  (68a)  does  not  affect  the  accuracy 
of  the  numerical  integration  with  respect  to  u. 


(71) 


-!l'  .»#" 
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where  u  ,  u  ,  and  u  are  unit  vectors  in  the  x,  y,  and  z  directions, 

-x  -y  —z 

respectively.  In  (71),  the  letter  t  stands  for  transmitter.  The 
origin  is  on  the  axis  of  the  body  of  revolution  and  is  in  the  vicinity 
of  the  body  of  revolution  although  not  necessarily  centered  Inside  it. 

In  (69)  and  (70),  the  E/ s  are  electric  fields  and  the  H’s  are  magnetic 
fields. 

Each  of  the  incident  plane  waves  (69)  and  (70)  comes  from  the 

direction  for  which  8  =  8fc  and  <j>  =  0.  Here,  8  and  <j>  are.  standard 

spherical  coordinates.  8  is  measured  from  the  positive  z  axis.  <j>  is 

the  angle  that  the  projection  onto  the  xy  plane  makes  with  the  positive 

t  t 

x  axis.  At  the  aspect  (8  =  8t»  <j>  =  0),  the  unit  vectors  u^  and  u^  re¬ 
duce  to  the  unit  vectors  in  the  8  and  $  directions,  respectively.  Be¬ 
cause  of  the  circular  symmetry  of  the  body  of  revolution,  no  generality 
is  lost  by  confining  to  the  xz  plane.  If  k^  were  rotated  through  an 
angle  4>t  from  <|>  =  0  to  <{>  =  <j>  ,  the  response  would  also  rotate  through 
the  same  angle  (j)^. 

In  view  of  (70a),  substitution  of  (69b)  for  H1  in  (11)  gives 

2tt 

I*f  dt  Tt(t)  |  d<|>  (Uj>*Ei<!>)e“jn(|>  (72a) 

0 

2tt 

~~  dtpPi(t)  d<|>  (^  \Ei<!>)e"jn<t>  (72b) 

1  0 

The  additional  superscript  8  on  the  left-hand  sides  of  (72)  indicates 
the  0  polarized  incident  plane  wave.  In  view  of  (69a),  substitution 
of  (70b)  for  H*  in  (11)  gives 

2tt 

~  |  dt  T±(t)  |  d<j>  (u^  •  Eie)e"jn4>  (73a) 

0 

2tt 

w~  dtpPi(t)  d<J>  (u^.- E19)e'jn<t>  (73b) 

1  0 
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The  third  superscript  on  the  left-hand  sides  of  (73)  indicates  the  <j> 
polarized  incident  plane  wave. 

Substitution  of  (9)  for  the  W's  and  (69a)  for  _E~  in  (7)  and  (8) 

of  [1]  reveals  that  the  elements  V*"®  and  considered  in  Section  IV 

ni  ni 

of  [1]  are  given  by 


t8  _  1 

ni  n 


dt  T±(t)  d<j>  (Uj.*  Ei0)e"jn(!> 


(74a) 


=  7ir  dtp?i(t)  d<^  (vii0)e 

1  r\ 


i9\  “jn<j> 


(74b) 


Substitution  of  (9)  for  the  W’s  and  (70a)  for  12  in  (7)  and  (8)  of  [1] 
gives 


,t<j>  _  1 


Vni  n  dtTi(t)  d*  >e 


,i<k  -jn<t> 


(75a) 


Vni  =  7^7  |  dt  P  Pi(t)  {  d<}> 

1  A 


(75b) 


Comparison  of  (72)  with  (75)  shows  that  if  T^(t)  and  (pP^(t))/p^  were 
interchanged  in  (75),  then 


it0  _  <M> 
ni  ni 


(76a) 


i<j>e  _  t(j> 

ni  “  ni 


(76b) 


Comparison  of  (73)  with  (74)  shows  that  if  T^(t)  and  (pP^(t))/p^  were 
interchanged  in  (74),  then 


i“*.  v*® 

ni  ni 


(76c) 


i#  =  _  te 

ni  ni 


(76d) 


Now,  T^(t)  is  responsible  for  both  the  asterisk  on  the  left-hand  sides 
of  (120)  and  (122)  of  [1]  and  the  factor 


25 


,  (-DP  2(t-t  ) 

In  + - 2_ 

2U  A 


) 


on  the  right-hand  sides  of  (120)  and  (122)  of  [1].  Also,  (pP  (t))/p 
appears  as  the  factor 

(t-t  )sin  v 
1  + - 2 - 2- 


in  (121)  and  (123)  of  [1].  Hence, 

n+l_,  r 2+1  (-l)p-i  2 (t-t  ) 


*ite  =  _  j 


ni 


ii<f>0  =  -  j\k 
np 


(1  + 


jkz  C08  0 

A - 


dt 


(77a) 


P+1 


(1  + 


(t-t  )sin  v 


5 - R)(Vl+  Jn-l)3lnVp  C 


jkz  cos  9 


dt 


*it<J>  =  j  TTk 
ni  2 


P+1 


(1  + 


(-1)1"1  2(t-t  ) 


+  Jn-l)oos  6t  6 


jkz  cos  0 


(77b) 


dt 


(77c) 


_.i<f>4>  ,n  . 
I  “  -  j  irk 
np 


P+1 


(t-t  )sin  v 

(1  + - ^ - ^)  [ j  sin  vp  cos  0t  ( -  Jn_,) 


-2  cos  v  sin  0,J  ]e 
p  t  n 


jkz  cos  0 


dt 


(77d) 


where  Jn»  p,  and  z  are  given  by  (115),  (130)  and  (13.1)  of  [1],  The 

asterisk  on  the  left-hand  sides  of  (77a)  and  (77c)  indicates  that  the 

right-hand  sides  of  (77a)  and  (77c)  are  not  all  of  and  but 

ni  ni 

only  the  contributions  due  to  the  region  of  integration  for  which 
tp  £  t  £  tp+i*  Both  (77a)  and  (77c)  are  valid  for  i  =  p  -1  and  i  =  p. 

Equation  (77)  can  be  rewritten  as 
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.1*4*1-  A 
♦ItS  3  rtA 
ni 


.  £  (F  , .  -  F  )  - 

4  n4*l , a  n*-l , a 


-  , xp-i.n-t-l  ,  . 

(-1)  j  kA 

- - - E.  (y  _  y  ) 

4  irKL,b  n-l,a 


C78a) 


jn 7r  kA  sin  v 

I1*6  -  -  - - ? - 2  [F  +  F 

np  2  -J-1 


A  sin  v 

.  *  -f  — ^ - E  (F  +  F  )1 

rrf-l,a  n-l,a  2p  v  n+l,b  n-l,b'J 


(78b) 


.  .  ,  jniT  kA  cos  8  (-1)P  ijnF  kA  cos  9 

fit* - E - £.  (.  +  F  )  +  - - - 2 - 1 

ni  4  n*HL,a  n-l,a  4  n+l,b  n-l,b 

(78c) 


ci^4>  _ 

np 


j’^irkA  sin  v  cos  8 


E - JL_ - E  r  _  p  )  -) — E-  -  — — E.(y  -f  )  1 

2  u  ttfl.a  n-l,a'  2p  v  nfl,b  n-l,b;j 


A  sin  v 
P 


A  sin  v 

+  jnirkA  cos  v  sin  0  [F  +  -2- — • — E  F  ,  ] 

J  p  p  t1-  na  2p  nb 

1? 


(78d) 


where  the  F's  are  given  by  (128)  and  (129)  of  [1] .  The  point  Gaussian 
quadrature  formulas  for  them  are 


“T  (nj 
Fma  '  l  Ai  Jm(k0i  sln  Ve 

X-=l 


jkz^cos  0fc 


"t  (nT)  (n.,i 

Fmt  =  J  h  xt  Jm(tflBin  Ve 


jkz^cos  @t 


m=n~l,n,n+l 


(79a) 


(79b) 


where  J  is  the  Bessel  function  of  the  first  kind  and  where  f3„  and  z„ 

(nT) 

are  given  by  (134)  and  (135)  of  [1].  The  abscissas  x0  and  weights 
(nT) 

are  tabulated  in  Appendix  A  of  [10]  for  several  values  of  n^,. 


m  ' 
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V.  NUMERICAL  RESULTS 


Computed  results  for  the  electric  currents  induced  by  an  axially 
incident  wave  on  a  sphere,  a  cone-sphere,  and  a  finite  cylinder  are  pre¬ 
sented  in  Figs.  1-14.  For  axial  incidence,  0  is  either  0®  or  180°  and 
the  only  non-zaro  excitation  vectors  for  the  6  polarized  plane  wave  (69) 
are 


fit 

X-1 

\ 

-1 

i 

It  is  evident  from  (12)  and  (13)  that 


y“ 

Ytt 

Yf 

-Yf* 

-l 

-1 

1 

1 

y*5 

Y# 

-Yf 

L-i 

-1_ 

Li 

1_ 

(80) 


(81) 


In  consequence  of  (80),  (81),  and  (10),  the  only  non-zero  column  vectors 

?C  and  are  given  by 
n  n  °  J 


“tt“ 

X-1 

a 

X1 

t* 

-1 

1 

(82) 


where  the  column  vector  on  the  right-hand  side  of  (82)  satisfies  (10) 
for  n  *  1. 

In  view  of  (3),  substitution  of  (82)  into  (2)  and  subsequent 
division  by  k  give 


J  „  t  _  i  P,(t) 

-j7l  ■ 2  ♦<!  rij  > +  2Vln  * (|  hi 

J 

The  |H|  written  instead  of  k  on  the  left-hand  side  of  (83)  is  the 


(83) 
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magnitude  of  the  incident  magnetic  field.  This  magnitude  is  indeed 
equal  to  k.  At  t  ®  fcp+i»  component  of  (83)  reduces  to 
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!e_ 


Ih1!  WKt^) 


cos  4>,  p=l,  2,  . .  .?-2 


(84) 


where  p(t  . ,)  is  the  value  of  p  at  t  =  t  ...  At  t  =  t  ,  the  <{>  component 

p+l  pr+-l  p 

of  (83)  reduces  to 


J,  2jlf  sin  $ 

- ,  P-  1,2,.  ..P-1 


kp„ 


(35) 


Here,  and  are,  respectively,  the  t  and  <j>  components  of  the 
electric  current  _J.  In  Figs.  1-14,  — —  in  the  <j>  =  0°  plane  is  plotted 

iri 

ijj 

with  squares  and  — J-  in  the  <j>  =  90°  plane  is  plotted  with  octagons. 

|H  I 

These  currents  are  plotted  versus  t/X  where  t  is  the  arc  length  along 
the  generating  curve  and  X  is  the  wavelength.  The  horizontal  axes  in 
Figs.  1-14  were  labeled  T/X  because  the  lower  case  letter  t  could  not 
be  drawn  by  the  plotter. 

Figures  1-4  display  computed  values  of  electric  current  induced 
on  a  conducting  sphere  of  radius  .2X  illuminated  by  a  plane  wave.  Figure  1 
shows  the  H-field  solution  of  [2],  Fig.  2  shows  the  present  H-field  solution, 
Fig.  3  shows  the  E-field  solution  of  (2],  and  Fig.  4  shows  the  E-field  solu¬ 
tion  of  [1].  The  squares  and  octagons  represent  computed  values  of 

|Jti  |JJ 

— r-  and  — ,  respectively.  The  solid  curves  represent  the  Mie  series 

Is1!  Ii1! 

solution  [13],  In  Figs.  1  and  3,  the  squares  and  octagons  are  horizontally 
located  at  the  peaks  of  the  triangles  used  [2].  In  Figs.  2  and  4,  the 
squares  are  placed  at  the  peaks  of  the  triangles  (5)  and  the  octagons  are 
placed  at  the  centers  of  the  pulses  (6) .  That  is  why  the  squares  and  the 
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Electric  current  on  a  conducting  sphere  of  radius  0.2a.  The 
squares  and  octagons  represent  the  E-field  solution  of  [2] . 
The  solid  curves  represent  the  Mie  series  solution. 
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Electric  current  on  a  conducting  sphere  of  radius  0.2X. 

The  squares  and  octagons  represent  the  E-field  solution 
of  [1].  The  solid  curves  represent  the  Mie  series  solution. 
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octagons  line  up  with  each  other  in  Figs.  1  and  3  but  are  staggered  in 
Figs.  2  and  4.  The  oscillations  of  tne  octagons  about  the  solid  curve 
in  Fig.  3  are  not  surprising  because  the  particular  size  of  the  sphere 
places  the  electric  current  error  in  Fig.  5  on  page  32  of  [2]  near  one 
of  its  peaks. 

Figures  5-10  display  computed  values  of  electric  current  induced 
on  a  cone-sphere  by  plane  waves  incident  from  both  axial  directions.  The 

uel  w 

squares  represent  — —  and  the  octagons  represent  — f- .  These  symbols 

Is  I  ll1! 

are  connected  by  straight  lines  in  order  to  improve  readability.  The 
currents  in  Figs.  5  and  6  agree  reasonably  well  with  those  in  Fig.  7  of 
[1],  Likewise,  the  currents  in  Figs.  8  and  9  agree  reasonably  well  with 
those  in  Fig.  8  of  [1].  The  jagged  nature  of  the  <j>  components  of  current 
in  Figs.  7  and  10  is  probably  spurious. 

Finally,  Figs.  11-14  display  computed  values  of  the  electric  cur¬ 
rents  induced  on  a  closed  cylinder  of  length  .5X  and  radius  .25X  by  an 
axially  incident  plane  wave. 

The  currents  in  Figs.  2,4,6,9,12,  and  14  were  calculated  with 


The  currents  in  the  rest  of  the  figures  were  calculated  with  n^  =  20 
in  [2].  The  IBM  System  370/155  under  WATFIV  took  29  seconds  cf  exe¬ 
cution  time  to  calculate  the  H-field  solution  in  Fig.  6.  The  sub¬ 
routine  YMAT  was  used  for  this  calculation.  The  joint  calculation  of 
the  H-field  solution  in  Fig.  6  and  the  E-field  solution  of  [1]  in 
Fig.  7  of  [1]  took  45  seconds  of  execution  time.  The  subroutine  YZ 
was  used  for  this  calculation.  The  joint  calculation  of  the  H-field 
solution  of  [2]  in  Fig.  5  and  th  E-field  solution  of  [2]  in  Fig.  7 
took  39  seconds  of  execution  time  using  the  subroutine  YZ  of  [14] . 
Hence,  the  present  H-field  solution  and  the  E-field  solution  of  [1]  are 
slightly  slower  than  the  H-field  and  E-field  solutions  in  [2], 
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Electric  current  induced  on  a  cone-sphere  by  an  axially 
incident  plane  wave.  Incidence  on  sphere,  H-field 
solution  of  [2]. 
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Electric  current  induced  on  a  cone-sphere  by  an  axially 
incident  plane  wave.  Incidence  on  sphere,  present  H-field 
solution. 
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Fig.  7.  Electric  current  induced  on  a  cone-sphere  by  an  axially 

incident  plane  wave.  Incidence  on  sphere,  E-field  solution 
of  [2]. 
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Fig.  9.  Electric  current  induced  on  a  cone-sphere  by  an  axially 
incident  plane  wave.  Incidence  on  tip,  present  H-field 
solution. 


incident  plane  wave.  Incidence  on  tip,  E-field  solution 
of  [2]. 
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The  results  of  calculation  of  the  present  H-field  solution  and 
E-field  solution  of  [1]  for  the  sphere,  cone-sphere  and  finite  cylinder 
examples  with 

\  "  *1 

were  so  close  to  the  results  in  Figs.  2,4,6,9,12,  and  14  that  distinc¬ 
tion  was  impossible. 
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PART  TWO 

COMPUTER  PROGRAMS 

I.  INTRODUCTION 

A  computer  program  which  implements  the  present  H-field  solution 
is  described  and  listed  in  Part  Two.  This  program  consists  of  the  sub¬ 
routine  YMAT,  the  function  BLOG,  the  subroutines  PLANE,  DECOMP,  and 
SOLVE,  and  the  main  program  for  the  H-field  solution.  The  subroutine 
YMAT  calculates  the  elements  (12)  of  the  moment  matrix  in  (10) .  The 
function  BLOG  is  called  by  YMAT.  The  subroutine  PLANE  uses  (78)  to 
calculate  the  excitation  vector  on  the  right-hand  side  of  (10)  for  a 
0  polarized  incident  plane  wave.  The  subroutine  PLANE  also  calculates 
the  excitation  vector  for  a  <j>  polarized  plane  wave,  but  this  vector  is 
not  used  in  the  main  program.  The  subroutines  DECOMP  and  SOLVE  solve 

the  matrix  equation  (10)  for  and  1^.  The  main  program  for  the  H-field 

n  n 

solution  obtains  the  electric  current  induced  on  a  conducting  body  of 
revolution  by  an  axially  incident  plane  wave.  The  subroutines  YMAT  and 
PLANE  are  designed  not  only  for  axial  incidence  but  also  for  oblique 
incidence.  For  axial  incidence,  n  =  +1  and  8^.  of  (71)  is  either  0°  or 
180°.  The  subroutines  YMAT  and  PLANE  admit  n  =  M1,M1+1, . . ,M2  where 
M2  _>  Ml  2l  0.  The  subroutine  PLANE  also  admits  arbitrary  values  of  0  . 

Formulas  on  pages  (28)  and  (29)  of  [2]  obviate  calculation  of  moment 
matrices  and  plane  wave  excitation  vectors  for  negative  values  of  n. 

A  computer  program  which  implements  both  the  present  H-field 
solution  and  the  E-field  solution  of  [1]  is  also  described  and  listed 
in  Part  Two.  This  program  consists  of  the  previously  mentioned  sub¬ 
programs  BLOG,  PLANE,  DECOMP,  SOLVE,  a  new  subroutine  YZ,  and  the  main 
program  for  both  the  H-field  and  E-field  solutions.  The  subroutine  YZ 
calculates  the  moment  matrices  in  (10)  of  the  present  report  and  in 
(6)  of  [1].  The  subroutine  YZ  calls  the  function  BLOG.  The  subroutine 
PLANE  calculates  the  excitation  vector  on  the  right-hand  side  of  (10) 
for  the  8  polarized  plane  wave.  PLANE  also  calculates  the  elements 
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and  -V^  of  (74)  for  use  on  the  right-hand  side  of  (6)  of  [1] .  In 
ni  ni 

addition,  PLANE  calculates  excitation  vectors  for  the  <j>  polarized  plane 
wave  but  these  vectors  are  not  used  in  the  main  program.  The  subroutines 
DECOMP  and  SOLVE  solve  the  matrix  equations  (10)  of  the  present  report 
and  (6)  of  [1].  The  main  program  for  both  the  H-field  and  E-field  solu¬ 
tions  obtains  both  the  present  H-field  solution  and  the  E-field  solution 
of  [1]  for  the  electric  current  induced  on  a  conducting  body  of  revolu¬ 
tion  immersed  in  an  axially  incident  plane  wave. 


II.  THE  SUBROUTINE  YMAT 


With  regard  to  (10),  the  subroutine  YMAT (Ml,  M2,  NP,  NPHI,  NT,  IN, 

EH,  ZH,  X,  A,  XT,  AT,  Y)  puts  in  Y  the  moment  matrices  Y  defined  by 

n 


Y  - 
n 


n  -  10,  M1+1....M2 


(86) 


Storage  in  Y  is  such  that  the  ith  column  of  Y  goes  from  Y((i-1)*N  + 

n 

(n-Hl)*N*N+l)  to  Y(i*N+(n-Ml)*N*N)  where 


N  »  2*NP  -  3 


(87) 


The  only  output  argument  of  YMAT  is  Y.  The  rest  of  the  arguments 

of  YMAT  are  input  arguments.  The  input  argument  IN  generalizes  YMAT 

far  use  with  problems  other  than  conducting  body  problems.  The  term 

(ttA  )/p  in  (22d)  and  the  quantities  in  (23)  are  multiplied  by  IN 
q  q  n 

in  the  subroutine  YMAT.  Thus,  the  argument  IN  should  be  1  for  the 
present  H-field  solution.  In  this  solution,  the  magnetic  field  is  evalu¬ 
ated  just  inside  S.  Magnetic  field  evaluation  just  outside  S  can  be 
obtained  by  setting  IN  ■  -1. 


Each  of  the  input  arguments  except  Ml,  M2,  and  IN  represents  a 
variable  in  Part  One  of  the  text.  The  correspondence  is  tabulated  as 


40 


Input 

Equation 

Argument 

Text 

Variable 

Number 

NP 

P 

3 

NPHI 

n4> 

36 

NT 

nt 

35 

RH 

kp(t^), 

3  =  1,2,. 

..P 

21 

ZH 

kz(t^) , 

(V 

x£ 

'  3  =  1,2,. 

..P 

21 

X 

£  =  1,2,. 

•  *n. 

<j> 

37 

/V 

A£ 

A 

£  =  1,2,. 

•  *nc|> 

36 

(nj 

XT 

xr 

*=> 

II 

h-1 

s* 

to 

“nt 

35 

(O 

AT 

A£-  ’ 

V  =  1,2,. 

••nt 

35 

Here,  p(t.)  and  z(t,)  are  the  values  of  p  and  z  at  t=t,  on  the  generating  curve 
3  j  3 

and  k  is  the  propagation  constant.  The  generating  curve  starts  at  t  =  t^ 

and  ends  at  t  =  t  .  The  input  argument  P  must  be  a  positive  integer  not 
less  than  3.  The  input  arguments  NPHI,  NT,  X,  A,  XT,  and  AT  are  Gaussian 
quadrature  data.  Lines  11  and  12  put  c^  and  c^  of  (31)  in  CT  and  CP, 
respectively.  The  subroutine  YMAT  calls  the  function  BLOG  which  is  de¬ 
scribed  and  listed  in  Section  III  of  Part  Two. 

Minimum  allocations  are  given  by 

COMPLEX  Y (M3*N*N) ,  GA(NPHI) ,  GB(NPHI) ,  GC(NPHI), 

G1A(M3) ,  G2A(M3) ,  G3A(M3) ,  GIB (M3),  G2B(M3), 

G3B (M3) ,  GIG (M3) ,  G2C(M3),  G3C(M3) 

DIMENSION  RH(NP) ,  ZH(NP),  X(NPHI) ,  A(NPHI), 

XT(NT) ,  AT (NT) ,  RS(NP-l),  ZS(NP-l),  D(NP-l), 

DR(NP-l) ,  DZ(NP-l),  C2(NPHI),  C3(NPHI), 

C4(M3*NPHI),  C5(M3*NPHI),  C6(M3*NPHI), 

R2 (NT )  ,  Z2(NT),  R7(MT),  Z7 (NT) 


where  N  is  given  by  (87)  and 


M3  -  M2  -  Ml  +1 


The  elements  of  Y  are  calculated  from  (23)  and  (24)  with  G 

n  2  dot 

evaluated  as  in  Section  III  of  Part  One.  DO  loop  11  puts  in  C2 
and  4sin^  (~  <{>^)  in  C3.  <p^  is  given  by  (37)  and  4sin^(^<i>^>  is  needed 
in  order  to  evaluate  (15)  at  <J>  =  <t>£  *  With  regard  to  (36),  inner  DO 


loop  29  puts 


2  1 

irA^  *  sin  (j  <j>^)  costn^)  in  C4, 


it 

2  cos  <}>£  cos(n(j)&)  in  C5, 


(n , ) 

2  Aj  ?  sin  <J>£  sin(n<{>£)  in  C6. 


The  index  JQ  of  DO  loop  15  obtains  q  in  (23)  and  (24)  .  With 
regard  to  (22d),  line  69  puts  (IN)*(ttA  )/p  in  P2.  DO  loop  12  puts 
kp1  and  kz'  of  (29)  in  R2  and  Z2,  respectively.  DO  loop  12  also 
accumulates  the  products  of  IN  with  (23a),  (23b),  and  (23c)  in  P1A, 
P1B,  and  PIC,  respectively.  The  index  IP  of  DO  loop  16  obtains  p 
in  (24).  Lines  106-112  put  kdQ  of  (32)  in  D6.  If  either  case  1  of 
(30)  or  case  3  of  (51)  is  obtained,  branch  statement  26  sends  execu¬ 
tion  to  statement  41.  If  either  case  2  of  (39)  or  case  4  of  (60)  is 
obtained,  execution  proceeds  to  line  114.  Lines  114-235  calculate 

G  for  case  2  and  case  4.  Lines  237-308  calculate  G  for  case  1 
ma  ma 

and  case  3.  Here,  G  is  defined  by  (25).  In  cases  2  and  4,  the 

ma 

integration  with  respect  to  t1  is  done  first.  In  cases  1  and  3,  the 
integration  with  respect  to  <j>  is  done  first.  In  all  cases. 


G  .  is  stored  in  GmA(n-Ml+l) 
ml 


G  „  is  stored  in  GmB(n-Ml-f-l) 
m2 

G  _  is  stored  in  GmC(n-Ml+l) 
ml 


m  -  1,2,3 
n  =  Ml,  Ml+1 , . .  .M2 


With  regard  to  (41),  DO  loop  33  puts 


vv 

in 

GA(K) 

H2(V 

in 

GB  (II) 

VV 

in 

GC(K) 

K  —  1^2, ••  •  n . 


where  <j>  is  given  by  (37).  Line  123  is  based  on  (48).  DO  loop  35 
accumulates  the  H's  of  (49)  in  H1A,  H2A,  and  H3A.  DO  loop  37  accumu¬ 
lates  in  H1A,  H2A,  and  H3A  the  sums  on  £’  in  (44).  If  kR  _<  .5,  then 
the  approximation 


1 


w 


2kR 


_  kR  (kR)2  (kR)4  j_  (kR)2 

8  v  18  720  ; ~  3  K  10 


+  (kH) 


280 


is  invoked  in  order  to  avoid  excessive  roundoff  error.  Lines  152  to  174 


put  the  I’s  of  (47)  in  Wl,  W2,  and  W3.  I.  and  I0  are  even  in  t  and  I„ 

l  j  o  z 


is  odd  in  tQ.  These  three  I's  are  first  calculated  with  t  replaced  by 

|tQ|  and  then  the  sign  of  is  changed  if  tQ  is  negative. 

2  -52  2-52 

In  line  161,  d  is  compared  to  10  r  .  If  d  <  10  r  ,  then 

pq  -  pq 

little  confidence  can  be  placed  in  the  calculated  value  of  d2  because 
2  2 

r  and  t  are  very  close  to  each  other  in  (46b) .  If 
pq  o 

2  -5  2 

d  <  10  r 

pq 


<±A 

2  q 


then  statement  52  stops  execution  because  an  accurate  value  of  1^  can 


not  be  obtained.  However,  if 

,2 


<  10‘5r2  , 
pq 


then  accurate  calculation  of  1^  may  be  possible. 


If  t 


2  \  is 


appreciably  greater  than  d  ,  then  the  approximation 


.2,2 

k  d  r 


|t  i  A 

o  2  q 


U  -?Aq 


'^2  I' 

2k 


<ig-lv2 


<ig+iv2 


can  be  invoked.  Line  163  puts  the  right-hand  side  of  (90)  in  W4. 
Lines  166-169  calculate  the  logarithm  term  in  (47a)  according  to 


[log(vfr)] 


(w?  +  r  ) 
108  [(„^-+  r,)1 


>  wx  >  0 


(w  +  r,)(-w1  +  r  ) 

log  [ - ~ - — ],  w  <  0 

d 


where 


(91) 


(92) 


This  logarithm  is  stored  in  W. 

Nested  DO  loops  45  and  46  calculate  (50)  for  n  =>  Ml  +  M-l  where 
M  is  the  index  of  DO  loop  45.  The  index  K  of  DO  loop  46  obtains  i  in 
(50).  With  regard  to  (68a),  lines  222-228  put 


in  D8.  DO  loop  67  adds  D8  to  G^.  DO  loop  67  also  sets  G ^ 

and  G^  equal  to  zero.  These  four  G's  are  set  equal  to  zero  when  p  =  q 


because  the  exact  values  of  the  coefficients  by  which  they  are  multiplied 
in  (24)  are  zero  then. 


The  index  L  of  DO  loop  13  obtains  V  in  (35).  With  regard  to  (36), 
DO  loop  17  puts  G(u,4>  )  in  GA(K) .  If  (59)  is  not  true,  line  258  sends 

K. 

execution  to  statement  51.  Lines  259-277  calculate  the  terms  in  (57)  which 
do  not  involve  G(u,  (J)^).  DO  loop  62  accumulates 


A  ,2 


L  3  3 

£-1  k  aj 


H  . 

—  in 


D6 


Line  275  puts  the  coefficient  of  n  in  (57c)  in  D8.  Line  276  puts  the  portion 

of  (57a)which  does  not  involve  G(u,<J>.)  in  D6.  Line  277  puts  in  D7  the 

%  2 

portion  of  (57b)  which  involves  neither  G(u,(J)^)  nor  (n  +1).  DO  loop  32 
accumulates  G  ,  G2>  and  of  (36)  in  H1A,  H2A,  and  H3A,  respectively. 

The  index  M  of  outer  DO  loop  30  obtains  n  =  M  4-  Ml-1.  If  (59)  is  true, 
lines  295-29/  add  to  H1A,  H2A,  and  H3A  the  terms  in  (57)  which  are  not 
present  in  (36).  Lines  298-306  do  the  sums  with  respect  to  in  (35). 


Lines  309-385  use  the  previously  calculated  G  to  obtain  the 
contributions  (23)  and  (24)  to  the  elements  of  the  moment  matrices  Y  . 

The  index  M  of  DO  loop  31  obtains  n  =  M  +  Ml-1.  With  regard  to  lines 
348-355,  the  subscripts  K1  to  K8  for  Y  are  the  same  as  the  subscripts 
K1  to  K8  for  Z  in  Table  2  on  page  50  of  [1],  The  variables  UA,  UG,  UB, 
UH,  and  UF  incremented  for  p  =  q  in  lines  338-342  are  intended  for  Y(K1), 
Y(K2),  Y(K3),  Y(K4),  and  Y(K8  +  MT) ,  respectively.  Y(K8  +  MT)  is  re¬ 
served  for  (Y^)  of  (24d) . 
n  pq 


45 


OO'C  LISTING  OF  THE  SUBROUTINE  YM  AT 

002C  THE  SUBROUTINE  YMAT  CALLS  THE  FUNCTION  BLOG 

003  SUBROUTINE  Y  MAT (  Ml  . M2. NP. NPH I* NT , tN.RH, ZH, X. A , XT, AT. Y) 

0  04  COMPLEX  Yt  1600)  .U.H1 A.H2 A.H3A.GA (48) , GB 1 48) , GC ( 4 8 ) . H 1 8, H28.  H38 

005  COMPLEX  HlC.H2C.H3C,UA.UB,UC.U0.UE,UF,GlAt 101.G2AI  10)  .G3AI10) 

006  COMPLEX  G1B< 10) .G2BI 10) .G3BC 10 ) *G1C< l 0 ) .G2CI tO) . G3CI 10).CHPLX 

007  COMPLEX  UG.UH 

0  06  DIMENSION  PH  (43)  .ZHI43).  X{  48  )  ,  A  (  48  )  .XT  ( l  0  )  .  AT  (  10  )  » RS  (  42)  .  ZS(  42  ) 

009  DIMENSION  D (42) . DR( *2)  .  D Z{ 42 ),  C2 14 8) , C3 1  48) , C4(  200 > .C5(  200) 

0  10  DIMENSION  C6  (200)  .R2I10)  .Z2  (  10  )  .  R7  (10  )  .  Z7  (1 0  ) 

Oil  CT=2. 

012  CP=.  I 

0  13  DO  10  l=2.NP 

014  12=1-1 

0  15  RS(I2)=.5*(RH(  I  >  ♦RH{  12)) 

016  ZS(t2)=.5*(ZHCI)  +  ZH(I2>) 

0  1 7  Dl=.£*(RH(I  )-RH(  12)) 

018  02=. 5*( ZH( I )— ZH ( 12) ) 

0  19  0(  I2)=S0RT(01*D1402*D2) 

020  DR  I  I  2  )  =0 1 

021  0Z(I2)=02 

022  10  CONTINUE 

023  M3=M2-Ml+l 

0  24  M4=Ml — l 

025  P  12=1  .570796 

026  PP=9. 869604 

027  DO  U  K=  t.NPHC 

3  28  PH=P(2*(X(K)+1 .) 

029  C2(K)=PH4PH 

030  SN=S(N( .5*PH) 

031  C3(K)=4.*SN*SN 

332  Al=PI2*A(K) 

0  33  -  D  4=  »  5*A  1  *C  3  (  K) 

034  D5=A1*C0S(PH) 

035  D6=A  1*SIN(PH) 

036  M5=K 

037  DO  29  M=l,M3 

038  PHM=  (M4  +  M)*PH 

039  A2=COS(  PHM  ) 

0  40  C4 ( MS  )=D4*A2 

041  C5( M 5) =D 5*A  2 

042  C6(M5)=D6*SIN(PHM) 

043  M5=M5HNPHl 

0  44  29  CONTINUE 

045  U  CONTINUE 
046  PNl=  .78539824IN 

047  PN2=8.*PNl 

0  48  U=  (0  .*  I  »  ) 

049  RP=NP-1 

0  50  MT  =  M  P—  l 

051  N=MT  + MP 

052  N2N=MT*N 

053  N2=N*N 

0  54  JN=—  1  -N 

055  DO  15  JQ=l «  NP 

056  KQ=  2 

357  IF(JO.EO.l)  KQ=t 

058  IF(JO.EQ.MP)  KQ  =  3 

059  R1=RS(JQ) 

060  Zl=ZS(JQ) 


061 

01=01  JO* 

0  62 

D2=DR(JQ) 

063 

D3*0Z<  JO) 

0  64 

04=02/81 

065 

05=0 t/R) 

0  66 

SV=D2/0t 

067 

CV=03/Ot 

0  68 

Pt=PN  l*Dl 

069 

P2=PN2*05 

0  70 

P3=2.*Dt 

071 

P4=2,*04 

072 

P5=D4*D4 

0  73 

P6=0l *01 

0  74 

P7=P6*0  l 

0  75 

T6=CT*D1 

0  76 

T62=T6+D  l 

0  77 

T62=T62  *  762 

078 

R6=CP*R l 

0  79 

F62=R6*R6 

080 

PtA=0  . 

081 

PI  B=0  • 

082 

PtC=0. 

0  83 

00  12  L=1.NT 

084 

D6=XT(L> 

0  85 

R2  (L)  =Rl  +02*06 

0  86 

Z2(L)=Zt +03*06 

0  87 

D7=Pl* AT (L)/R2(L> 

088 

D8=  1  .-06 

089 

09=D8*07 

090 

P1A=P  IA +08*09 

0  91 

06  =1  .+06 

0  92 

P 1 B=P IB  +06*09 

093 

PI C=P 1 C+C6 *06*07 

0  94 

12 

CONTINUE 

0  95 

DO  16  t  P  =  t  »  MP 

096 

R3=RS(1P» 

0  97 

Z3=ZS(tP) 

0  98 

R4=R  l— R3 

0  99 

Z4=Z1-Z3 

too 

DO  40  L= 1  *  NT 

101 

D7=R2(LI-R3 

l  02 

D8=Z2(U-Z3 

103 

R71L  )=R3*R2(L> 

t  04 

Z7( U =D7*07+08*D8 

105 

40 

CONT  1NUE 

l  06 

PH  =R4*  SV+  Z4*  CV 

107 

A  l  =  A8S ( PH) 

l  08 

A2=ABS( R4*CV— Z4*  SV) 

l  09 

06  =  A2 

t  10 

tFlAl  .t_E.Dl )  GO  TO  26 

t  U 

D6= A l— D l 

t  12 

D6=SQRT(D6*D6+ A2*  A2) 

l  13 

26 

IF!  1P.NE.JQ  .AN0.1R6.GT.D6.0R  .T6.LE.06)  )  GO  TO  41 

l  14 

Z5=R4*R4+Z4*Z4 

1  15 

RS=R3*R1 

1  16 

PHM=.5*R3*SV 

l  17 

DO  33  K*  UNPHI 

t  18 

Al =C3  <  K ) 

l  19 

RR  =  Z  5+R5*A 1 

120 

HI A=0. 

47 


121 

H2A=0 . 

122 

H3A=0. 

t  23 

1F1RR.LT.T62)  GO  70  34 

1  24 

DO  35  L= 1 .  NT 

125 

W=Z7(L)+R7(L)*A1 

125 

R=SQRT( W ) 

127 

SN=-S(N(R) 

1  28 

CS  =  COS(  R ) 

1  29 

H 1 B= AT ( L  )/ 1 R*W  ) *CMPL  X( CS-R*SN, SN+R4CS) 

130 

HI  A=Hl S+Hl A 

131 

H2B=XT(L  )  *H  IB 

132 

H2A=H28+H2A 

t  33 

H3A  =  XT{ L  >*H2B+H3A 

t  34 

35 

CONT CNUE 

1  35 

GO  TO  36 

1  36 

34 

00  37  L=l .NT 

l  37 

W=Z7(L)+R7(L)*A l 

l  38 

R=SQRT (W  ) 

139 

IF(R  ,GT .  .5)  GO  TO  14 

t  40 

CS=R*  (W* 1 .  69444 44E-2-W*. 17361 1  IE-3)- .125) 

141 

SN  =  Vrf*( . 3233333E— l-W* . 1 190476E-2)-. 3333333 

142 

Hie=AT(L)*CNPLX(CS.SN> 

143 

GO  TO  43 

1  44 

14 

SN=-S 1N1R) 

1  45 

C S=C0 S(  R  ) 

1  46 

HI  3=  AT(L)/R*((CMPLX  <  CS— R  *SN  «  SN+R  *C5)—  l  .)  /W-  .5) 

l  47 

43 

H  l  A=H  IB  +  H  l  A 

148 

H2B=XT(L)*H1B 

t  49 

H2A=H2B+H2A 

1  50 

H3A=XT(L)*H2B+H3A 

151 

37 

CONTI  NOE 

1  52 

At=PH+PHM*At 

l  53 

A2=ABS(  At) 

154 

R=RR-A2*A2 

l  55 

D6=A2-D l 

l  56 

D7=A2+Dl 

1  57 

D62=D6*D6 

l  58 

D72=C7*D7 

t  59 

D  8  =  SORT (062+R) 

l  60 

D9=SCRT  ( 072+R) 

l  61 

(F ( R-(RP*l , E-5) )  52.52.53 

1  62 

52 

CF1D6.LT.0.)  STOP 

l  63 

te4=. 5/062-. 5/072 

l  64 

GO  TO  54 

l  65 

53 

W4  =  ( D 7/0 5—0 6/08 ) /R 

l  66 

54 

IF1D6.GE.0.)  GO  TO  38 

1  67 

W-AL0G1 (D7+D9)*(-06+C8)/R) 

1  68 

GO  TO  35 

1  69 

38 

#=ALGGC (07+09)/ (D6+D3 ) ) 

170 

39 

W 1  =  1  W4+  .5**1  )/D l 

171 

W5= A2/D1 

l  72 

W2  =  (  .5*1  09-08 >- l ./D9+1./D8) /P6- W5*  W 1 

l  73 

W3=(  .25*  (07*09-06*08  5 +  W-R* (W4+ .25*VO  >/P7-W5* ( 2 . *W 2+W 5*W l  ) 

l  74 

(F ( A  l  .LT  .0.  )  W2»-W2 

l  75 

HI  A=WI+HIA 

l  76 

H2A=  *2+H2A 

l  77 

H3  A=VI3  +  H3A 

t  78 

26 

GA(K )=HIA 

l  79 

GS(K)=H2  A 

l  80 

GC1K  )=H3A 

tat 

33 

CONTINUE 

t  82 

Kt=0 

133 

00  45  M=t ,M3 

t  34 

81A=0  . 

135 

H2A=0. 

l  36 

H3A=0  . 

tar 

Ht8=0. 

tas 

828=0 » 

t  39 

H36=0  . 

l  <50 

H  tC  =  0  • 

1 9  t 

H2C=0 . 

192 

H3C=0. 

t  93 

DO  46  K=t,NPHt 

194 

K l =K  t  + 1 

t  95 

D6=C4(  Kt  ) 

t  96 

D7=C5( Kl > 

1 97 

D8=C6  (Kt  ) 

1  98 

UA=G A { K) 

199 

Ue=GE(K) 

2  00 

UC=GC( K) 

20  t 

Ht  A=D6*UA+Hl A 

202 

H2A=0  7*U  A4-H2A 

203 

83  A=D8*U  A +  83 A 

2  04 

H  IB=D6*U8+H 1 0 

20  5 

82 8=07  4UB+H28 

206 

H3B=0a*U8+H3B 

207 

81 C=D6*UC+81C 

2  08 

H2C=D7*UC+H2C 

209 

83C=08*UC+H3C 

2  1 0 

46 

CONTINUE 

211 

Gl  At  M  )=H  t  A 

2  1 2 

G2A( N)=H2A 

213 

G3A( M }=H3A 

2  14 

GIBt  N) =H t  B 

215 

G2B(M  )=8£B 

2  16 

G3B( N)=H3B 

2  17 

GIC(M)=HIC 

210 

G2C( H)=H2C 

219 

G3C(M)=83C 

2  20 

45 

CONTINUE 

2  2 1 

tF(  tP.NE  .JQ)  GO  TO  47 

222 

At  =05*05 

223 

08=0  . 

224 

DO  63  K=l,NPHt 

225 

O8=O0*A(  K  )/SQRT(  C2(  K  )• 

2  26 

63 

CONTINUE 

227 

A2  =  3»l4t  593/05 

2  28 

D8=tetOGt A2  J-Pt2*D3)/ 

2  29 

00  67  M=l.M3 

?  30 

Gl A(Mt=08+Gl A(M) 

231 

G2A(M)=0. 

232 

G2E( M)=0  . 

233 

G2C(M)  =  0. 

2  34 

G3  A  (  M  t  =0  t 

2  35 

67 

CONTINUE 

236 

GO  TO  47 

237 

41 

DO  25  M=t.M3 

2  38 

Gl  A (  M )  =0  • 

239 

G2A( M  >  =0 • 

240 

G3  A  (  N  )  =0  . 
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\ 


I 


1 

I 

j 


241  GtB<M}=0. 

242  G2B(M)=0. 

243  G3B( M  )  =  0  • 

244  GlCCM)=0. 

245  G2CC  M )=0  < 

246  G3C(M)=0. 

247  25  CONTINUE 

248  DO  13  L=l.NT 

249  RS=R7(L) 

250  Z5  =  Z7  { L ) 

251  DO  17  K=t,NPH{ 

252  W=Z54R5*C3{tO 

253  R~SQRT(Vf) 

254  SN=— S  IN ( P) 

255  CS=C0S{R) 

256  GAIK  )=CMP1_X(CS-P*SN.SN4R*CS)/{W*R) 

257  17  CONTINUE 

258  tFlR62.LE.Z5)  GO  TO  St 

259  06=0. 

260  07=0  . 

26t  OO  62  K=  t . NPHt 

262  W2=C2(K) 

263  W=l./(Z54R5*W2t 

264  W  l=AIK)*SQRT<WJ 

265  D6=064Wt*W2** 

266  D7=D74W  V*(  .S4-W*<  1.4.  125* W*R5*W2*W2) ) 

267  62  CONTINUE 

268  W  t  =R 5/Z5 

269  N2=PP*W1 

270  W=S0RT(W2) 

271  *3=1, 4W2 

272  R=SORT { W  3) 

273  W4=S0RT(R5> 

274  -W 5=Al_0G(  W4-R  ) 

275  D8=— P 12*  C6— ( #/R— W5}/ (R5**4 ) 

276  06= . 5*08 

277  D7  =  < (W/R*(Wl-( ,t25+.  1666667**2 )/ W 3 )  +  . I 25*W 5)/R 5+ . 5*W 5) ZWA-P l 2*07 

278  51  Al=AT(L) 

279  A2=XT(L)*A1 

280  A3=X  Ti L  > *A 2 

281  K l =0 

282  00  30  M=t.M3 

283  *=M>M4 

284  H  t  A=0 . 

285  N2A=0. 

286  H3A=0. 

287  00  32  K  =  t . NPHI 

288  Kl=Kt4t 

289  HI E=GA( K ) 

200  H  IA=C4( Kl»*HtB4HlA 

29 1  H2A=C5(KtJ*HlB4H2A 

292  H3A=C6(  Kt)*HtB  +  H3A 

293  32  CONTINUE 

294  IF C R62.LE.Z5)  GO  TO  44 

295  HI A=D6*M l A 

296  H2A  =  0  7—  I  W*  W4  l  .  )  *064H2A 

297  H3  A=**08  4H3A 

298  44  GIA(  N)=Al*HlA4Gl  AIMI 

299  G2A(M)=Al*H2A*G2AIM) 

300  G3AOO=At*H3A*G3A(M) 


a 


| 

5'fl 

1 
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]  301 

G1 8(M)  =  A2*H1A4-G1B{M) 

361 

tFUP.EQ.MP)  GO  TO  2* 

1  302 

G28{  M  >  =  A  2*H  2A4-G  2B<  M) 

3  62 

21 

Y (  K4) =UH 

|  30  3 

G3B(M>  =  A2*H3A4-638(M> 

363 

Y  (  K8  >=UE 

'  304 

GIC( M)  =  A3*HtA4-61C<M> 

364 

GO  TO  22 

1  305 

G2C(M  }  =  A3*H2A4-G2CtM» 

365 

19 

Y(K5  )=Y  { i(5  )  ♦UC— UD 

l  3  06 

G3C  (  M )  =  A  3*H3A4-G3C  <  M> 

3  66 

IFdP.EO.U  GO  TO  23 

307 

30  CONTINUE 

367 

Y {  Kl  )=Y(KU+UA 

303 

J 

13  CONTINUE 

3  68 

Y1  K7)=Y«  K7)-rUE 

309 

47  A2=0 { IP } 

3  69 

»F(tP.EQ«MP>  GO  TO  22 

310 

A3=.S*A2 

3  70 

23 

Y<  K2)=Y(K2)+UG 

311 

W t=A 3 *(R4 *03-24 *D2f 

371 

Y(K8)=UE 

312 

W2=— A3*  R3*D3 

3  72 

GO  TC  22 

]  313 

A3=DZ ( IP  1 

373 

20 

Y(K5)=Y(K5)+UC-U0 

,  314 

D6=0R  ( t  P ) 

3  74 

Y£ K6)=UC4UD 

i  3  15 

D7=Z4  *06 

375 

1F(  tP.EQ.l)  GO  TO  24 

j  316 

09=0  3*06 

3  76 

Y(Kl »=Y( Kt )4UA 

j.  317 

HtC=(Rl *09— 02*( R  3*A  3+D7 J )*U 

377 

YC X3)=Y( K3)4U8 

j  3  18 

H3C=A2*0 1*U 

378 

Y1K7)=Y(K7)4UE 

j  319 

H2C= Z 4*H3C 

379 

1F(  IP  ,EQ  .MPJ  GO  TO  22 

;  3  20 

i 

H3C=03*H3C 

3  80 

24 

Y  (  K2  )  =Y  (  K2  )  4UG 

1  321 

W3=P3*(R4*A3-07> 

381 

Y(K4)=UH 

3  22 

V»4=P3*(  02*  A3— 09  ) 

382 

YCK8)=UE 

3  23 

W5=P3*Rl *A3 

383 

22 

Y( K84MT  5  =UF 

!  3  24 

JM  =  JN 

384 

JM=jm+N2 

3  25 

DO  3  1  M=  1,  M3 

385 

31 

CONTINUE 

3  26 

H2A=G2A<  M) 

386 

16 

CONTINUE 

3  27 

H l A=  G l A { M) 

387 

JN=JN4N 

,  3  28 

H  2  8=  62  81  M) 

3  88 

15 

CONTINUE 

3  29 

H1B=61B( M> 

3  89 

RETURN 

330 

UC=Wt  *H2 A+W2*H1  A 

390 

ENO 

331 

UB=Wl*H2B+W2*HtB 

332 

UF=W3*{H2A4-04*H28HW4*(H2B«-D4*G2C(M)  >4-W 5*( H  IA4-P4*H1B*P5*S ICC  M)  ) 

333 

UA=UC—  us 

3  34 

U0--UC4U8 

3  35 

UG=UA 

3  36 

UH=UB 

337  tF<  IP  «NE  .JQ  )  GO  TO  48 

3  38  UA=P  l  A+U  4  i 

339  UG=P  184-UG 

340  UB=Pte*UP 

3  41  UH=P IC4UH 

342  UF=P24UF 

343  48  H3A=G3A(N> 

344  H 38=G  39  ( M}  ; 

345  UC=HIC*H3A 

3  46  UD=HtC*H3B 

347  UF=H2C*(H3A404*H3B)4H3C*  ( H3B4D 4*G3C( M ) ) 

3  48  K!=IP4JM 

3  49  K2=K141 

350  K3=!<  l-t-N 

3  5 1  K4=K24N  , 

352  K5=K2+MT 

353  K6=K4  4MT 

354  KT=K34N2N 

355  K8=K  44N2N  j 

356  GO  TO  (1  8»20.19)  . KQ 

357  18  Y(K6)=UC4UD  ;i 

358  IF(tP.EO.l)  GO  TC  2.  j 

359  Y(K3)  =  Y(K3»4-UB  j 

360  Y(  K?1  *Y<  K7»  4-UE  1 
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Ill .  THE  FUNCTION  BLOG 


/  2 

The  function  BLOG(x)  calculates  log(x  +\J1  +  x  )  for  x  >  0. 
The  method  of  calculation  is  described  on  page  56  of  [1] . 


00  l  C  LtSTING  OF  THE  FUNCTION  BLOG 

002  FUNCTION  BLOG<  X ) 

003  t  F  (  X  •  GT  «  .t  )  GO  TO  t 

004  X2=X*X 

005  0LOG=  (  (  •  0754X2—  .  1 666667  )  *X  ?.+ 1  .  )  *X 

006  RETURN 

007  1  BCOG=ALOG(X+SQRT(l .+X*X) > 

008  RETURN 

009  END 
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IV.  THE  SUBROUTINE  PLANE 


The  subroutine  PLANE (Ml,  M2,  NF,  NP,  IN,  NT,  RH,  ZH,  XT,  AT, 
THR,  RE,  R)  obtains  the  plane  wave  excitation  vectors  (7 2)  —  (75) .  If 
IN  ^  2,  PLANE  stores  the  plane  wave  excitation  vectors  (72)  and  (73) 
for  the  present  H-field  formulation  in  RE.  If  IN  ^  1,  PLANE  stores 

Vtj,  -V^®,  -V^,  and  of  (74)  and  (75)  in  R.  The  minus  signs 
ni  nl  nl  ni 

attached  to  and  transform  the  V's  into  plane  wave  measurement 
ni  ni  r 

vectors.  Plane  wave  scattering  calculations  are  not  done  in  this 
report,  but  if  they  were  done,  plane  wave  measurement  vectors  would  be 
required.  I^6  is  put  in  RE(i+(n-Ml)*2*N+(K-l)*(M2-Ml+l)*2*N)  for 


N  =  2*NP-3 
i  =  1,2, ...NP-2 
n  =  Ml,  Ml+1,.  ..M2 
K  =  1,2, ...NF 

Ii(^  is  displaced  forward  NP-2  locations  from  1^®  in  RE.  For  I^f^, 
ni  ni  ni 

i  =  1,2,...NP-1.  1^^  is  displaced  forward  N  locations  from  1^®, 

ni  ni 

I^f^  is  displaced  forward  N  locations  from  I^f®.  The  angle  of  inci- 
Hj.  ni 

dence  of  (71)  is  THR(K)  radians  and  K  =  1,2,.,.NF.  The  storage 

arrangement  of  V^,  -V*^,  -V*^,  and  in  R  is  the  same  as  that  of 
ni  ni  ni  ni 

the  I’s  in  RE.  The  arguments  RE  and  R  of  PLANE  are  output  arguments. 
The  rest  of  the  arguments  of  PLANE  are  input  arguments.  Some  of  the 
input  arguments  obtain  variables  in  Part  One  of  the  text.  The  cor¬ 
respondence  is  tabulated  as 
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Input 

Argument 

Text 

Variable 

Equation 

Number 

NP 

P 

9 

NT 

79 

RH 

kp(t^), 

j  =  1,2, ...P 

21 

ZH 

kz(t  ), 
(nT) 

j  =  1,2, ...P 

21 

XT 

\nT)' 

£  —  1,2,... nT 

79 

AT 

A  1 
£ 

£  =  l,2,...nT 

79 

Here,  k  is  the  propagation  constant.  Also,  p(t..)  anc*  z(t^)  are 

values  of  p  and  z  at  the  point  t  =  t.  on  the  generating  curve.  This  curve 

starts  at  t  =  t.  and  ends  at  t  =  t  .  It  is  assumed  that  P  is  not  less 

1  r 

than  3.  The  input  arguments  NT,  XT,  and  AT  are  Gaussian  quadrature 
data. 


Minimum  allocations  are  given  by 

COMPLEX  RE(2*N*(M2-M1+1)*NF) ,  R(2*N*(M2~M1+1)*NF) , 
FA(M2+3),  FB(M2+3) 

DIMENSION  RH(NP),  ZH(NP) ,  XT (NT) ,  AT (NT) ,  THR(NF) , 
CS(NF) ,  SN(NF) ,  R2 (NT) ,  Z2(NT) 


The  I's  and  V's  are  calculated  from  (78)  of  the  present  report 
and  (124)-(127)  of  [1].  The  index  IP  of  DO  loop  12  obtains  p  in  these 
equations.  With  regard  to  (79),  DO  loop  13  puts  and  kz^  in  R2 

and  Z2,  respectively.  The  index  K  of  DO  loop  14  obtains  the  Kth  angle 


of  incidence.  Nested  DO  loops  15  and  25  accumulate  S*FW  „  and  S*F„  ^  , 

M-2,a  M-2,b 

in  FA(M)  and  FB(M),  respectively.  Here,  F^  ^  a  anc^  2  b  are  8^ven  by 
(79)  and  S  is  the  normalizing  constant  for  the  Bessel  functions.  Accord¬ 
ing  to  (9.1.46)  on  page  361  of  [15]  , 


iTi  v  nrf  fin  rtftw 
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S  =  j“(x)  +  2j“(x)  +  2J^(x)  +  ...  (93) 

where  the  superscript  u  means  unnormalized,  Jm(x)  are  the  cylindrical 
Bessel  functions  of  the  first  kind  and 


x  =  kf^  sin  (94) 

The  logic  inside  DO  loop  15  is  the  same  as  that  inside  '.he  DO  loop  15 
which  appears  in  the  listing  of  the  subroutine  PLANE  on  pages  61  and  62 
of  [1].  If  Ml  =  0,  lines  88  and  89  use  the  formulas 

(95) 

to  store  F  -  and  F  in  FA(1)  and  FB(1),  respectively. 

”"-L  >  B-  i-  >  D 

The  index  M  of  DO  loop  27  obtains  n  =  M-2.  In  DO  loop  27, 


I  ^  is  added  to  RE(Kl)  for  i  =  p-1 

*iffl 

1^^  is  put  in  RE(K2)  for  i  -  p 
I^Q  is  put  in  RE(K3) 

I^^^is  added  to  RE(K4)  for  i  =  p-1 


In^tl^is  put  in  RE(K5)  for  i  =  p 

1^^  is  put  in  RE(K6) 
np 

The  above  I's  are  given  by  (78).  The  V's  are  given  by  (124) — (127 )  of  [1]. 
Location  of  these  V's  in  R  is  similar  to  location  of  the  I's  in  RE. 
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00 t  LISTtNC  OF  THE  SUBROUTINE  PLANE 

002  SUBROUTINE  PL AN E ( Ml . M2 , NF. NP« I N, NT , RH, ZH, XT . AT , THR. RE.R ) 

003  COMPLEX  RE ( 240) «R ( 240)  .U.Ul.FA(lO) ,FB(  l 0 ) ,UA, UB . F l A , F t B . F2A ♦ F2B 

004  COMPLEX  U2.U3.U4.U5.CMPLX 

005  0  IMENSION  RH(43)  .2H(  43 ) , XT < 10 > , AT ( 10)  ,  THRO  >  .  CS<  3  ).  SN<3  )  »R2  (  10  ) 

006  DIMENSION  Z2  (  10  )  .  BJ  I  50  ) 

007  MP=NP— 1 

008  MT=MP-I 

009  N=MT+MP 

0  10  N2  -2  *N 

Otl  00  l  l  K=t,NF 

012  X=THR(K) 

013  C  S ( K) =C0  SIX) 

014  SN  (  K  )=S  I  N  (X  ) 

015  U  CONTtNUE 
0  16  U=  (0  .tl  .  ) 

017  Ul=3.  l4tSS3*U**Ml 

0  18  M3  =  Mt+l 

019  M4=M2+3 

020  IF(Ml.EQ.O)  M3=  2 

0 2  t  M5  =  Ml  +  2 

0  22  M6=M2  *2 

023  00  12  IP=l.MP 

024  K2=IP 

025  t=IP+l 

026  DR=.5*(RH(t )-RH< IP)) 

027  DZ=. 5*1 ZHI I )-ZH ( tP) ) 

028  0  l=S QRT  I  0R4DR+-DZ  *0Z  ) 

029  Rt =.25* ( RH( l )+RH ( IP) ) 

030  Zl=.5*<ZH(l)+ZH< tP)) 

031  02=.5*0R 

032  DR=D2/R1 

033  DO  13  L  =  1  *  NT 

034  R2<L1=RI+D2*XT(L> 

035  Z2(L) =Zl +OZ*XT(L ) 

0  36  13  CONTINUE 

037  W 1 =— . 5  *  D l 

038  W2=—  2  «*D 2 

039  DO  14  K=l,NF 

040  CC=CS(K) 

041  SS=SN(K) 

042  03=02*CC 

043  04=-DZ*SS 

044  O5=-0t*CC 

045  W 3-— 2  .*0  3 

046  W4  =— 2 •*  04 

047  W  5=— • 5*0  5 

048  00  23  M=M3tM4 

049  FA  (  M  }=0  • 

050  F8(M)=0. 

051  23  CONTINUE 

052  00  15  L= l ♦ NT 

053  X=  SS4R2I L ) 

054  I F (X . GT  • »5E“7 )  GO  TO  19 

055  00  20  M=M3»M4 

056  P J ( M )  =0  • 

057  20  CONTINUE 

058  BJ(2)=l, 

0  59  S=l. 

060  GO  TO  18 
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061 

19 

M=2.8*X+  14.-2. /X 

121 

U2=D  3*F l A4U5*F  A  <H) 

062 

IFIX.LT..5)  M=l  1  .8+ALOGl 0 IX ) 

l  22 

U3=D3*F 1 6+US*FB(M> 

063 

IF1M.LT.M4)  M=M4 

123 

U4=D  27F  2 A 

064 

B  J 1 M  )  =0  . 

124 

U  5=D24F2  B 

0  65 

JM=M- l 

t  25 

IF(iP.EQ.l)  GO  TO 

066 

8J(  JK)  =  t  . 

1  26 

R ( X l  )  =  R {K 1 ) 4U2— U3 

067 

DO  16  J=4,M 

127 

R( K4>  =  R< K4) +U4-U5 

0  68 

J2=JM 

128 

tF(tP.EQ.MP)  GO  TO 

069 

JM= J  M— l 

t  29 

30 

Rl K2)=U2+U3 

0  70 

Jl =JM-t 

1  30 

R ( K5 )=U  44U5 

071 

BJ< JM)=J l/X*BJ (J2)— 8JIJM+2) 

1  31 

29 

K2=K2+N2 

0  72 

16 

CONTINUE 

132 

UA=UB 

073 

s=o. 

1  33 

27 

CONTINUE 

074 

IF1M.LE.4)  GO  TO  24 

1  34 

14 

CO  NT  tNUE 

0  75 

00  17  J=4.M.2 

1  35 

12 

CONTINUE 

0  76 

S=SF6J( J) 

136 

RETURN 

0  77 

17 

CONTINUE 

137 

END 

078 

24 

S=ej(2)+2.*S 

0  79 

18 

ARG=Z2(C)*CC 

0  80 

U A  =  AT (L ) /S* CMPLX (COS(ARG). SIN(  ARG) ) 

081 

UB=XT{ L ) *UA 

0  82 

DO  25  M=M3.M4 

0  83 

F  A(  M  )  =8  J  1  M  )T  UA4-F  A(  M) 

084 

FBIM  )  =  BJ(M>*UB  +  FB(M) 

0  85 

25 

CONTINUE 

0  86 

15 

CONTINUE 

0  87 

IF1MUNE.0)  GO  TO  26 

0  88 

F  A  (  l  )=-FA(3l 

0  89 

FBI U --FBI3) 

0  90 

26 

UA=Ul 

091 

DO  27  M-N5.M6 

092 

M7=M-l 

093  M  8=M ♦ J 

094  F2A=UA*(FA(M8)+FAIM7)> 

095  F2B=UA*(FB( M8)+F8<M7)) 

096  UB=U*UA 

097  FIA=UB*( FA( M8)-FA1M7)) 

098  FI  e=UB*(F6(M8)-FB(M7)) 

0  99  K  l  =K  2- 1 

100  K3=K2+MT 

tOt  K4=Kl+N 

102  K5=K2+N 

103  K6=K3+N 

104  IF(  tN.EQ.2)  GO  TO  28 

105  RE(K3)=W2*(F2A+DP*F2E) 

106  RE1K6  )=W3*(FtA  +  DR*FlS)+W4<>UA*(KA(MM-DR*FB<MH 

107  U2=W  t*F  t  A 

1 08  U3=WI*FIB 

t  09  U4=W5*F2A 

l  10  U5=W5*F2  S 

t  t  t  IF (l P.E  0 .  t )  GO  TO  2 1 

t  t  2  REIKI  )=RE(K  l»«-U2-U3 

t  t  3  RE<K4)=RE<K4)*U4-U5 

t  14  IF(IP.EQ.MP)  GO  TO  22 

ttS  21  RE(K2)=U2-M»3 
t  16  RE  (K 5  )-U  •.  f 

t  17  22  t F ( IN.EQ.l )  ,  0  TO  29 

118  28  R(K3  )  =  054(F2A  '>R 4F2B  ? 

t t9  R(K6)=D1*(F1Av'-lR*F1B> 

120  U  5=0  4*U A 
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V.  THE  SUBROUTINES  DECOMP  AND  SOLVE 

The  subroutines  DECOMP(N,  IPS,  UL)  and  SOLVE(N,  IPS,  UL,  B,  X) 
solve  a  system  of  N  linear  equations  in  N  unknowns.  The  input  to 
DECOMP  consists  of  N  and  the  N  by  N  matrix  of  coefficients  on  the  left- 
hand  side  of  the  matrix  equation  stored  by  columns  in  UL.  The  output 
from  DECOMP  is  IPS  and  UL.  This  output  is  fed  into  SOLVE.  The  rest  of 
the  input  to  SOLVE  consists  of  N  and  the  column  of  coefficients  on  the 
right-hand  side  of  the  matrix  equation  stored  in  B.  SOLVE  puts  the 
solution  to  the  matrix  equation  in  X. 

Minimum  allocations  are  given  by 

COMPLEX  UL(N*N) 

DIMENSION  SCL(N) ,  IPS(N) 
in  DECOMP  and  by 

COMPLEX  UL(N*N) ,  B(N) ,  X(N) 

DIMENSION  IPS(N) 

in  SOLVE. 

More  detail  concerning  DECOMP  and  SOLVE  is  on  pages  46-49 


of  [16] . 


ootc 

Lt STING  OF  THE  SUBROUTINES  OECOMP  AND 

SOLVE 

002 

SUBROUTINE  CECOHP< N. CPS, UL 1 

003 

COMPLEX  UL( 16001 .PC VQT,EM 

004 

0 1 MENS  ION  S CL (40  )* IPS(40 ) 

005 

00  5  t=l,N 

006 

1PS( I )*t 

007 

RN=0. 

008 

Jt=( 

009 

DO  2  J=l , N 

0  10 

ULM=A8S(  PEALCULC  J 1 > > l+A&SC A IM AG ( UL ( J  111  >■ 

0  1 t 

Jl=Jt+N 

0  12 

IF(RN-ULM)  1.2,2 

013 

1 

RN=ULM 

0  14 

2 

CONT  INUE 

0  15 

SC  L(  t)  =  t./RN 

0  16 

5 

CONT  INUE 

017 

NMt=N-l 

0  18 

K2=0 

019 

DO  17  K=l,NNl 

0  20 

BtG=0. 

021 

DO  lt  t  *K»N 

0  22 

IP=  IPSC I  1 

0  23 

IPK=tP+K2 

0  24 

StZE=  CA8SCREALCULC  IPX  1  ))+ABS(AlMAGC ULC IPKl )> )*SCLC  Cr) 

0  25 

IF(SIZE-etG)  ll ,11.10 

026 

10 

8  IG=S  IZE 

0  27 

CPV=t 

028 

U 

CONT INUE 

0  29 

IF(tPV-K)  14.15,14 

0  30 

14 

J=  IPSCK) 

031 

tPSCK)=IPS( IPV) 

032 

IPS!  CPV>=J 

0  33 

15 

KPP=  t  PS ( K  >  +K2 

034 

PlVOT=UL(KPP) 

035 

KPl=K+l 

036 

DO  16  t=KPt,N 

0  37 

KP=KPP 

038 

IP=  IPS(  !  14-K2 

039 

EM=-UL( tPl/PXVOT 

0  40 

ie 

UL  (  I P  )=—EM 

041 

OO  16  J=KPt,N 

061 

DO  l  J=l,lMt 

0  42 

IP=IP+N 

062 

SUM=$UM+ULCIP)*XCJ» 

0  43 

KP=KP+N 

0  63 

l 

IP=IP4N 

0  44 

ULC IP)  =  UL( I P )+E W4ULC  KP) 

064 

2 

X(n=B(IPB)-SUM 

045 

16 

CONT  INUE 

0  65 

K2=N* (N— 1 ) 

046 

K2=K2+N 

066 

l P=t  PSC  N)4K2 

047 

17 

CONTINUE 

0  67 

XCN)=X(N)/UL(tP> 

048 

RETURN 

066 

DO  4  I BACK=2 . N 

0  49 

END 

069 

l=NP  l-  I  BACK 

050 

SUBROUTINE  SOLVE ( N, I  PS ,UL. E, X ) 

0  70 

K2=K2-N 

051 

COMPLEX  UL ( 1600 1 .B(40),X{ 401. SUM 

071 

IPt=  tPSC  t)+K2 

0  52 

DIMENSION  I  PSC40) 

072 

tPt=t+l 

0  53 

NP1=N+1 

073 

SUM=0. 

0  54 

1P=( PSC  1 1 

0  74 

!P=tPl 

0  55 

X( l )  =  8  C IP) 

075 

DO  3  J=  t PI , N 

0  56 

DO  2  1=2, N 

076 

IP=IP4N 

057 

tP= IPS (  I  1 

077 

3 

SUM=SUM+UL  C IP)*XC J) 

058 

I  PB  =  I P 

0  78 

4 

xc  n  =  cx(  t >— sum )/ul c  tPii 

0  59 

CM  1=  t—  l 

079 

RETURN 

060 

SUM=0. 

080 

ENO 
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VI.  THE  MAIN  PROGRAM  FOR  THE  H-FIELD  SOLUTION  : 

The  main  program  for  the  H-field  solution  calculates  the  H-field 
solution  for  the  electric  current  induced  on  a  conducting  body  of  revolu¬ 
tion  immersed  in  an  incident  plane  wave.  Input  data  are  read  according 
to 

READ (1, 15)  NT,  NPHI 

15  FORMAT (213) 

READ ( 1,10) (XT (K),  K  -  1,  NT) 

READ (1,10) (AT(K) ,  K  -  1,  NT) 

10  FORMAT (5E14. 7) 

READ(1, 10) (X(K) ,  K  =  1,  NPHI) 

READ(1,10) (A(K) ,  K  =  1,  NPHI) 

READ(1,16)  NP,  BK,  THR(l) 

16  FORMAT (13,  2E14.7) 

READ(1,18) (RH(I) ,  1=1,  NP) 

READ(1,18) (ZH(I) ,  1=1,  NP) 

18  FORMAT (10F8. 4) 

The  input  data  obtain  variables  in  Part  One  of  the  text.  The  input  data 
are  tabulated  versus  text  variables  in  the  following  chart. 
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Input 

Data 

Text  Variable 

Equation 

Number 

NT 

nt 

35 

NPHI 

<V 

xr 

36 

XT 

9 

£'  •=  1,2,. . .n^ 

35 

AT 

/V 

V 

9 

£*  =  1,2, . . .n^ 

35 

X 

(V 

xj£ 

» 

l  =  1,2,...!^ 

37 

A 

A(V 

» 

£  =  1,2, ...t^ 

36 

NP 

P 

3 

BK 

k 

69 

THR(l) 

9t 

71 

RH 

p(tj> 

9 

j  =  1,2, ...P 

21 

ZH 

2(V 

.  9 

j  =  1,2, ...P 

21 

(nt)  (n  )  (n  )  (n  ) 

The  Gaussian  quadrature  data  x^,  ,  A^,  ,  x^  Y  ,  and  are  given  in 

Appendix  A  of  [10].  In  the  main  program  for  the  H-field  solution,  n^  of 
(79)  is  set  equal  to  n  .  NP  controls  the  order  of  the  moment  matrix  Y 

^  X 

of  (86).  This  order  is  (2*NP-3) .  THR  is  in  radians.  THR  should  be 
either  0  or  tt  depending  on  the  direction  from  which  the  incident  wave 
approaches.  The  subscript  on  THR  is  for  compatibility  with  the  subroutine 


PLANE.  P(tj)  and  z(t^)  are  t^ie  values  of  p  and  z  at  the  point  t  *  tj  on 
the  generating  curve  of  the  body  of  revolution.  This  curve  starts  at 


t  =  t^  and  ends  at  t  =  tp.  p  is  the  distance  from  the  axis  of  the  body 
of  revolution  and  z  is  the  coordinate  measured  along  this  axis.  RH  and 


ZH  are  in  meters. 


The  main  program  for  the  H  field  solution  calls  the  subroutines 
YMAT,  PLANE,  DECOMP,  and  SOLVE.  The  function  subprogram  BLOG  is  also 


needed  because  it  is  called  by  the  subroutine  YMAT. 


Minimum  allocations  are  given  by 

COMPLEX  Y(N*N) ,  R(2*N) ,  C(N) 

DIMENSION  XT (NT) ,  AT (NT) ,  X(NPHI) ,  A(NPHI) , 
RH(NP),  ZH(NP) ,  IPS (N) 

where 

N  =  2*NP-3 


The  t  and  $  components  of  the  electric  current  are  calculated  from 
(84)  and  (85).  The  coefficients  I^p  and  I^  in  these  equations  are  the 
pth  elements  of  the  vectors  1^  and  1^  which  satisfy  the  n  =  1  equation 
in  (10) .  DO  loop  28  prepares  RH  and  ZH  for  use  in  the  subroutines  YMAT 
and  PLANE  by  multiplying  them  by  k.  With  regard  to  (10)  line  41  puts 
Y^  of  (86)  in  Y.  Line  46  calculates  IPS  and  changes  Y.  Line  47  puts 
the  excitation  vectors  1^  and  1^  for  the  0  polarized  incident  plane 
wave  (69)  in  R.  Line  47  also  stores  the  excitation  vectors  for  the  <}> 
polarized  incident  plane  wave  (70)  further  on  in  R  but  these  vectors 
are  not  used  in  the  main  program.  In  line  50,  the  output  IPS  and  Y  from 
the  subroutine  DECOMP  is  fed  along  with  N  and  R  into  the  subroutine 
SOLVE,  SOLVE  puts  and  in  C. 


The  pth  line  printed  out  under  the  heading  JT  contains  the  real 
part,  the  imaginary  part,  and  the  magnitude  of  the  normalized  t  component 
(Jt)/|l±l  electric  current  at  <J>  =  0°  and  t  =  t^^.  The  pth  line 
printed  out  under  the  heading  JP  contains  the  real  part,  the  imaginary 
part,  and  the  magnitude  of  the  normalized  tj;  component  (d(j))/|H1|  of  electric 
current  at  <j)  =  90°  and  t  =  t^.  According  to  (69b), 

lii1!  -  -  4  •  [s\-o 


where  [H^]z_q  is  the  incident  magnetic  field  at  z  =  0.  The  sample  input 
and  output  data  are  for  the  sphere  example  of  Fig.  2.  The  input  array 
ZH  was  constructed  so  that  z  =  0  at  the  center  of  the  sphere. 
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OOtC  LISTING  OF  THE  MAIN  PROGRAM  FOR  THE  H-FIELD  SOLUTION 

002C  THE  SUBPROGRAMS  YMAT , BLO G, PL ANE, DECOMP ♦  AND  SOLVE  ARE  NEEDED 

003/ /PGM  JOB  ! XXXX.XXXX, l • 1) ,*MAUTZ,JCE* .REG! ON=20OK 

004//  EXEC  WATptV 

005 / /GO*  SYS l N  DO  * 


006$ job  MAUTZ.T  IME=5,  PAGES*  60 

007  COMPLEX  Y ( 1 600)  * R (240) «C! 4  0 )  »U  »C  t 

008  0  IMENS  ION  XT(10).AT(  10).  X!  48).A(  '  8)  .RH C 43) ,ZH( 43) . THR< 3)  «tPSt40) 

009  READ ( 1.15)  NT.NPHC 

0  10  15  FCRMAT(2(3> 

Ott  WRITEI3.30)  NT. NPH! 

0  12  30  FORMAT ( *  NT  NPH t •/ IX,  13. I S) 

013  READ! t, 10) <XT(K) ,K=l .NT) 

014  READ! 1. 10)!AT!K),K=1.NT) 

015  10  FORMAT! 5E 14  «7) 

0  1 6  WRtTE(3»  11)(XT(K),  <=  1  .NT  ) 

0  t  7  WR  t  TE 13  »  12  )  ! AT <  K)  ,K=l .NT ) 

018  It  FORMAT!*  XT*/! IX.5E14.7) ) 

019  12  FORMAT!*  AT » / ! 1 X ,5El 4 .7 ) ) 

020  READ! t, 10) ! X!K) .K=t.NPHt) 

021  READ! I , 10 ) ! A!K) »K=l ,KPHl ) 

022  WRITE!3,  13)(X(K)  »K«  1  ,NPH C ) 

023  WRITE!3.t4)!A!K),K=i,NPHt) 

024  13  FORMAT!  •  X « /! I X . 5E X 4. 7) ) 

025  14  FORMAT!*  A* / ! 1 X ,5E1 4 .7 ) ) 

026  READ!  I,  l  6)  NP.8K.THRU) 

027  16  FORMAT! T3, 2614. 7) 

028  WR  IT E( 3 .  17 )  NP.BK.THRll) 

029  17  FORMAT!*  NP* ,6X . * BK * . t2X « 'T HR ' / IX. t 3. 2E 14.7) 

030  READ!  1,  l€)(RH!  I )  .1=1 ,NP) 

931  READ!!  ,1  8)  !  ZH{  I).  1=1  ,NP) 

032  16  FORMAT!  10F8.4) 

033  WR(TE(3»19)!RH!  D.I-l.NP) 

034  WR  ITE!  3.20MZH!  !).t  =  t,NP) 

035  19  FORMAT!*  RH*/!IX,10F8.4)  ). 

036  20  FORMAT! »  ZH*/! 1X.10F6.4) ) 

037  00  23  J=l,NP 

038  RH(  J  )=BK*RH(  J) 

039  ZH  !  J  )  =  RKvZH!  J) 

040  28  CONTINUE 

041  CALL  YMAT! l.l.NP. NPH l,  NT .  l.RH.ZH.X.A, XT  » AT«  Y) 

042  MT=NP— 2 

043  N=2*  MT+l 


044 
945 
046 
047 
048 
049 
050 
0  51 
0  52 
0  53 
054 
9  55 
0  56 
0  57 
053 
0  59 
0  60 


WP  ITE! 3,29 )( Y! J) ,J=l ,N) 

29  FORMAT!*  Y* /11X.661 l  .4)) 

CALL  DEC  OMP !N,!PS,Y) 

CALL  PL ANE!t .1  *  1 .NO. l.NT.RH, ZH, XT, AT ,T HR ,R, R ) 
WR  ITE!  3,  23)!R!  J)  ,J=1  ,N) 

23  FORMAT!'  R'/! IX.6E1 1  .4)) 

CALL  SOLVE! N, IPS, Y.P.C) 

U=(0  .  .1  .) 

WR  1  TE!  3 , 2  l  ) 

2 *  FORMAT! •  REAL  JT  IM AG  JT  MAG  JT*) 

DO  24  J=l,MT 

Ct  =2  ,/RH!  J+l )*C<  J) 

C£=CA5S! Cl ) 

WP!TE<3,25)  CI.C2 
25  FORMAT!  1X.3EU.4) 

24  CONTINUE 

WR 1 TE( 3, 26) 


)1j*? 
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06 1  26  FORM  AT  (  •  REAL  UP  I M  AG  JP  MAG  jo • ) 

062  MP=NP-t 

0  63  00  27  J=l,MP 

0  64  Cl=4./(RH< J  )+RH< J+t>)*U*C<  J+MT) 

065  C2=CA8S(Ct) 

066  WR ( TE( 3  .  25)  CI.C2 

067  27  CONTINUE 

0  68  STOP 

069  END 

SDATA 
2  20 

-0.57735 03E+00  0 . 5773503E+00 
0 .10000  08  E+0 l  0.  lOOOOOOE+Ol 

-  0.993t286E+00-0.96397l  9E+00-0 .9122344E+00— 0  .839t  170E+00-0  .7463319E+00 
-0 • 6 36 05 37 E  +00-0.5108670E+00-0  .  3 73 70 61 E+ 00-0 . 2 277859E+ 00-0 . 7652652 E-Ol 

0 • 765  26  52E— 0 1  0. 2277 859E+00  0 .373706 l E+00  0  .5 t 03670E+00  0 .6360537E+00 
0 .7463319E+00  0 .8391  l  70E+00  0  .9  122344E+ 00  0.  96397  19E+00  0 .9931 286E+00 
0. l76t40tE-0l  0. 40 60143 E- 01  0  .626 72 05 E-0 1  0 .8327674E-0  l  0 . 10 1930 t E+00 
O  .  I  181945E+00  0. 1316886E+00  0  .  I  420961E+00  0. 1 491 730E+00  0.15275345+00 
0.  1527534E+00  0.  149 1 73OE+0O  0 . I 42096t E+00  0  . I  3 16886E+00  0.118 1945E+00 
0 . 1019301E+00  0 . 8327674E- 0 1  0 .6267205E-01  0. 4060 1 43E-0 1  0 . 1 76 1 401 E-Ol 
16  0. 1256637E+01  0. 3 l 4 l 593E+0 t 

0.0000  0.2079  0.4067  0.5878  0.7431  0.8660  0.9511  0.9945  0.9945  0.9SU 

0.8660  0.7431  0.5878  0.4067  0.2079  0.0000 

-1.0000  -0.9781  -0.9135  -0.8090  -0.6691  -0.5000  -0.3090  -0.1045  0.1045  0.3090 

0.5000  0.6691  0.8090  0.9t3S  0.9871  1.0000 

i  STOP 
/* 

// 

PRINTED  OUTPUT 
NT  NPH  I 
2  20 
XT 

-0. 5773503E+00  0 .5773S03E+00 
A  T 

0. 1000000E+01  0.1000000E+01 
X 

-  0.993  l  286 E+00— 0  .96397  1  9E+00— 0  .9122344  .'*•-<)  .8391  1  70E  +  OO—O  »  7463319E+00 

-  0.6360537E  +00—0 . 5  108670E  +  00— 0 .3737061,.  .2 27785 9E+00-0  .7652652E-01 

0. 7652652E-0  t  0 .2277859E+00  0  .373706  1  £+Cj  ,  5  108670E  +  00  0 . 6360537E+00 
0. 746331 9E +00  0 . 839 1  1  70E  +  00  0 . 9  12  2344E+00  0 . 963  97 1  9E+00  0  .993 *  286E+00 

A 

0 , 1 76  140 1 E-0 l  0 .4060  143E-0  1  0  .6267208E- 01  0. 8327675E- 0 S  0 . 1 0 1 9301 E+00 
0. t  181945E  +  00  0.  131 6886E  +  00  0 . 1 42 096 1  E+00  0 . 1 4 91 730E+00  0 . 1527534E+00 
0 . 1527534E+00  0  . I  49 l 7 30E+00  0 . 1  4209 6 1 E  +  00  0 . I  3 1 6886E+ 00  0.1181 945E  +  00 
0. 1019301E  +  00  0.8327 675 E-Ol  0  .6267206 E-0 l  0 .4060 143E-0  l  0 . 176  1 40 1 E-Ol 
NP  BK  THR 

16  0. I 256637E+01  0. 3 l 4 l 593E+0 l 
RH 

0.0000  0.2079  0.4067  0.5878  0.7431  0.8660  0.9511  0.9945  0.9945  0.951 t 

0.8660  0.7431  0.5878  0.4067  0.2079  0.0000 

ZH 

-1.0000  -0.9781  -0.9135  -0.8090  -0.6691  -0.5000  -0.3090  -0.1045  0.1045  0.3090 

0.5000  0.6691  0.8090  0.9135  0.9871  t.0000 

Y 

0. 24  86E  +  01-0.233  8E-0  2  0  .42  It E  +  00-0 .49 I9E-1  2  0  .6397E-0 1-0.8888E-02 
0  .559QE-0  1-0 .  138  2E-P  1  0 . 4 S e8E- 0  t-0.  19  235-01  0 .44  55E-01-0  .2463 E-0 1 
0.394  IE-01-0 .2964E-0 l  0  .3  432E-0  t-0 .339BE-01  0 .2937E-0 1-0.3752E-01 
0  .24  75E-0  t-0  .  <0  2  35-0  l  0 . 2065S-  0  l-  0 .  42  175-01  0.  1  72  85-0  l -0 . 4348  E-0  1 
0. 1475E-01-0.4458E-0 t  0 . t 306E-0 1-0 . 4490 E-0 l  0.9666E-03  0.9243E-01 
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-0  .9622E-03-0 

.  5803E-0  1-0 

4589E-0 2-0. 6791E— 01-0. 951 2E-92-0.5939E-01 

-0. 151  7E-01-0 

.5297E-0 1-0 

2105E-0 l— 0 .4  766 E— 01  — 0  .2668E-0 1-0 . 4258E-01 

-0  .3164E-0 1-0 

.3741E-0  t- 

0 

3578E— 0  l- 0 . 32 26E— 01  —  0. 3903E— 0 t-0s2732E-0l 

-0.4137E-01-0 

•  22  79E— 0 1  - 

■0 

4297E-0 1-0 .1891 E-0t-0 .4402E-01-0. 1588E-0t 

-0  .4484E-0 1—0 

.  1362E-01-0 

4470c— 01— 0» 12  44E— 0 1 

K 

-0> 2834E+00-0 

•  7  74  eE4  00— 0 

3446E 400-0.74916400-0. 43 36E400-0 . 6972 E400 

-Of  5523E400-0 

•6097E400-0 

6671 E400-0 .4300  C400-0  .76095*00-0. 3088E400 

-0  .81 36E+00— 0 

.  t067E400-0 

8 1 36E400  0. t067E400-0.7609E400  O.3O83E40O 

-0.6671E+00  0 

.4800E  400  -0 

5523E400  0 .60 97 E400 -0.4386E400  0.6972E400 

-0  » 34S4E  +  00  0 

. 7550E400— 0 

2791E400  0.78 10E 400  0 .7758E400-0. 2679E4O0 

0. 7267E400-0 

.2959E400 

0 

6295E 400-0. 3378 E4i)0  0 . 489 1E400-0.  368AE400 

0  •32-36E40O— 0 

.3603E+00 

0 

1626E 400-0. 2921 E400  0 . 4374S-0 1-0 . 1652E4O0 

O.OOOOE+OC  0 

.0000E400 

0 

4374E-01  0.16S2E400  0.1626E400  0.2931E4OO 

0 • 3236E 400  0 

•  36G3E40  0 

0 

489  tE  400  0»  3684E4Q0  0.6295E400  0.33T8E400 

0.728 tE+00  0 

.2922E400 

0 

7778E400  0.26  20E400 

REAL  JT 

1 MAG  JT 

MAG  JT 

-  0  » 798 9 E +00—0 

. 1959E+0  l 

0 

21 15E40  l 

-  0.9480E400-0 

. t930E40l 

0 

21S1E40  t 

-0  .1 1  62E40  1—0 

. 1836E+0  1 

0 

21 73E4C1 

-0.  t406E  +  0  l—  0 

. 1658E40  t 

0 

21 75E401 

-  0  » 1630E+0 l— 0 

. 1372E+0  t 

0 

2130E401 

-0.  1  764E  +  0  t— 0 

.  S802£v00 

0 

2018E401 

--0 .17S4E401-0 

.5127E400 

0 

1827E40  t 

-  0. 1S69E  +  01-0 

.  242  8E-01 

0 

t  569E401 

-0 .1218E+0  l  0 

.4228E400 

0 

1289E401 

-  0. 7475E400  0 

.  7772E400 

0 

1078E401 

-0.2294E400  0 

.  !0  16E40  1 

0 

1042E401 

0.259 2E 400  0 

.  1 146E+0  l 

0 

l 1 74E40 1 

0.6S69E400  0 

. 1210E+0  l 

0 

1377E401 

0.9016E400  0 

.  1 190E40  1 

0 

1493E401 

REAL  JP 

IM4G  JP 

MAG  JP 

0.773  IE4-00  0 

.  1878E4  01 

0 

203 1 E40 1 

0. 8555 E 400  0 

. 1782E40  l 

0 

1976E40 l 

0.9955E400  0 

.  1602E40  1 

0 

1888E401 

0 . 1 1 48E40 1  0 

.  1333E40  1 

0 

1759E401 

0.1240E401  0 

.9997E400 

0 

1593E401 

0.1229E401  0 

«  C565E+0  0 

0 

13S4E 401 

0.U06E401  0 

.3761E400 

0 

1  168E401 

0.9065E400  0 

.2217E40C 

0 

9332E400 

0.7010E400  0 

.  2213E400 

0 

7351E400 

0 . 5589E400  0 

. 2543E400 

0 

6618E400 

0.5201E400  0 

.56316400 

0 

7665E400 

0.5789E400  0 

.7805E400 

0 

972  IE  400 

0.6903E400  0 

•  553  8E400 

0 

l  177E4©  1 

0.8141E400  0 

.  1078E40  l 

0 

1351E401 

0.8958E400  0 

.U4BE40  t 

0 

1459E40  i 

VII.  THE  SUBROUTINE  YZ 


With  regard  to  (10),  the  subroutine  YZ(M1,  M2,  NP,  NPHI,  NT,  IN, 

RH,  ZH,  X,  A,  XT,  AT,  Y,  Z)  puts  in  Y  the  moment  matrices  Y  defined 

n 

by  (86).  With  regard  to  (6)  of  [1],  the  subroutine  YZ  puts  in  Z  the 

moment  matrices  Z  defined  by 
n  J 


n  =  Ml,  Ml+1,  ...  M2 


(96) 


The  Y  obtained  by  YZ  is  exactly  the  same  as  the  Y  obtained  by  the  sub¬ 


routine  YMAT  which  is  described  and  listed  in  Section  II  of  Part  Two. 


The  storage  arrangement  of  the  elements  of  Z^  in  Z  is  the  same  as  the 
arrangement  of  the  elements  of  Y^  in  Y.  The  Z  obtained  by  YZ  is  nearly 
the  same  as  the  Z  obtained  by  ZMAT  of  [1].  The  greatest  difference  in 
calculation  is  that,  in  accord  with  (68a),  the  second  sum  with  respect 
to  £'  in  (97)  of  [1]  is  replaced  by  the  exact  value  of  the  correspond¬ 
ing  integral  with  respect  to  t'.  The  purpose  of  the  subroutine  YZ  is 
to  obtain  Y  and  Z  more  efficiently  than  by  means  of  YMAT  of  the  present 
report  and  ZMAT  of  [1] .  Increased  efficiency  is  possible  because  several 
calculations  done  in  YMAT  of  the  present  report  are  repeated  in  ZMAT  of  [ 1] . 


The  only  output  arguments  of  YZ  are  Y  and  Z.  The  rest  of  the  argu¬ 
ments  of  YZ  are  input  arguments.  The  input  arguments  of  YZ  have  the  same 
names  and  meanings  as  the  input  arguments  of  YMAT.  Lines  13  and  14  put 
cfc  and  c^  of  (31)  in  CT  and  CP,  respectively.  The  subroutine  YZ  calls  the 
function  BLOG  which  is  described  and  listed  in  Section  III  of  Part  Two. 

Minimum  allocations  are  given  by 

COMPLEX  Y(M3*N*N),  Z(M3*N*N),  GA(NPHI) ,  GB(NPHI) , 

GC(NPHI) ,  GD(NPHI) ,  GE(NPHI) ,  G1A(M3) ,  G2A(M3) , 

G3A(M3) ,  GIB (M3) ,  G2B(M3),  G3B(M3),  G1C(M3), 

G2C(M3),  G3C(M3),  G4A(M3) ,  G5A(M3) ,  G6A(M3) 

G4B(M3) ,  G5B(M3),  G6B(M3) 
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DIMENSION  RH(NP),  ZH(NP) ,  XCNPHI),  A(NPHI) ,  XT (NT) 

AT (NT) ,  RS(NP-l),  ZS(NP-l) ,  D(NP-l),  DR(NP-l) 

DZ(NP-l),  DM(NP-l),  Cl(NPHI),  C2(NPHI),  C3(NPHI) 
C4(M3*NPHI),  C5(M3*NPHI),  C6(M3*NPHI),  R2(NT) , 

Z2 (NT) ,  R7 (NT) ,  Z7 (NT) 

Here,  N  and  M3  are  given  by  (87)  and  (88),  respectively. 

The  elements  of  Y  are  calculated  from  (23)  and  (24)  with  G 

n  ma 

evaluated  as  in  Section  III  of  Part  One.  The  elements  of  Z  are  calcu- 

2  n 

lated  from  (48)-(51)  of  [1].  DO  loop  11  puts  <K  in  Cl,  <j>?  in  C2,  and 

2  1  ^  2  1  * 

4  sin  (y  4^)  in  C3.  is  given  by  (37)  and  4  sin  (-^  tp^)  is  needed  in 

order  to  evaluate  (15)  at  $  =  <}>£.  With  regard  to  (36),  inner  DO  loop  29 


(n.)  2  1 

irA^  9  sin  (y  <j)^)  cos(ncj>£) 


in  C4, 


1  * V 


cos  4>2  cos(n<j)^) 


\  a£  ^  sin  cj)£  sin(n(J>^) 


in  C5, 


in  C6. 


The  index  JQ  of  DO  loop  15  obtains  q  in  (23)  and  (24)  of  the  present 

report  and  in  (48)-(51)  of  [1].  With  regard  to  (22d) ,  line  75  puts 

(IN)*(irA^)/p^  in  P2.  DO  loop  12  puts  kp'  and  kz'  of  (29)  in  R2  and  Z2, 

respectively.  DO  loop  12  also  accumulates  the  products  of  IN  with 

(23a),  (23b),  and  (23c)  in  P1A,  P1B,  and  PIC,  respectively.  The  index 

IP  of  DO  loop  16  obtains  p  in  (24)  of  the  present  report  and  in  (48)-(51) 

of  [1].  Lines  114-120  put  kd  of  (32)  in  D6>  If  either  case  1  of  (30) 

o 

or  case  3  of  (51)  is  obtained,  branch  statement  26  sends  execution  to 

statement  41.  If  either  case  2  of  (39)  or  case  4  of  (60)  is  obtained, 

execution  proceeds  to  line  122.  Lines  122-288  calculate  G  for  case  2 
r  ma 

and  case  4.  Lines  290-388  calculate  G  ^  for  case  1  and  case  3.  In  cases 

ma 

2  and  4,  the  integration  with  respect  to  t'  is  done  first.  In  cases  1 
and  3,  the  integration  with  respect  to  <p  is  done  first.  In  all  cases, 
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G  .  is  stored  in  GmA(n  -  Ml+1) 
mi 


G^  is  stored  in  GmB(n  -  Ml+1) 


G  „  is  stored  in  GmC(n  -  Ml+1) 
mo 


G  is  stored  in  GmA(n  -  Ml+1) 
ma 


G  ,  is  stored  in  GmB(n  -  Ml+1) 
mb 


m  =  1,2,3 
n  =  Ml,  Ml+1,... M2 


(97) 


m  =  4,5,6 
n  =  Ml,  Ml+1,... M2 


(98) 


The  G's  in  (97)  are  defined  by  (25).  The  G's  in  (98)  are  defined  by 
(56)- (57)  of  [1]. 

With  regard  to  (41),  DO  loop  33  puts 


H1((J>k)  in  GA(K) 


H2(<j>K)  in  GB(K) 


K  1,2, .. . n , 


H3  (<J>k)  in  GC(K) 


DO  loop  33  also  evaluates  G  and  G  of  (71)  —  (7 2)  of  [1]  at  (j)  =  <j>  and  puts 

them  in  GD(K)  and  GE(K),  respectively.  The  branch  statement  in  line  133 

is  based  on  (48).  DO  loop  35  uses  (49)  to  calculate  H  (4>  )  of  (41).  The 

Ot  K. 

index  L  of  DO  loop  35  obtains  £'  in  (49).  DO  loop  35  uses  the  Gaussian 
quadrature  formulas 


Ga  =  l 

£’=1  *  kR 


Gk  =  l  A 


?  /V  (V  e"jkR 


Jb  g,L .  V  kR 

x  =1 


to  calculate  G  and  G,  of  (71) — (72)  of  [1],  In  (99), 

3  b 


(99a) 

(99b) 


* ,  ,  j  A'  V  jVhf*  '■ '  -  Vw 
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R  =  /(p'-pp) 2  +  (z’-zp)2  +  4ppp'sin2(-|  <J>K) 


where 


p’  =  P  + 


Vr  slnvc 


z  1  =  z  + 


<V 

x. ,  A  cos  v 
V  q  c 


(100) 


(101a 


(lOlt 


DO  loop  35  accumulates  H-  (<J>  )  ,  H0(<j>  ),  H_  (<j)  ),  G  ,  and  G  in  H1A,  H2A, 

H3A,  H4A,  and  H5A,  respectively.  DO  loop  37  accumulates  the  sum  on 
V  in  (44)  for  a  =  1,  2,  and  3  in  H1A,  H2A,  and  H3A,  respectively.  DO 
loop  37  accumulates  and  G  ^  of  (79)-(80)  of  [1]  in  H4A  and  H5A, 
respectively.  If  kR  _<  .5  where  R  is  given  by  (100),  then  lines  154-159 
use  (89)  and  ~ 

e~jkR-l  „  kR,_  1  ,  (kR)2  _  (RR)\  +1  /  ,  .  (kR)2  _  (kR)4  (kR)6 

kR  ~kR(-  2  +  24  -  720  )+j(-l+  6  -  12Q  +  5040  )  (102) 


in  order  to  put 


(\)  p  ! 

V  (G  -  -^3  -  2^)  inHIB 
(kR)  ^ 


(nt)  ,e-jkR-l, 


AV  kR~>  in  H4B 


If  kR  >  .5,  lines  161-165  calculate  H1B  and  H4B  directly  from  (103). 

Lines  173-199  put  1^  of  (47)  for  a  =  1,2,  and  3  in  Wl,  W2,  and  W3, 

respectively.  Lines  173-199  also  put  G  .  and  G  of  (87)-(88)  of  [1] 

a  z  bz 

in  W4  and  W5,  respectively.  1^,  1^,  and  G^2  are  even  in  tQ.  I,,  and 

G,_  are  odd  in  t  .  I, ,  I„ ,  I„ ,  G  and  G,  „  are  calculated  with  t 

b2  o  12  3  a/  b2  o 

replaced  by  |tQ|  and  then  the  signs  of  I2  and  G^2  are  changed  if  tQ  is 
negative.  Line  173  puts  kt_  in  Al. 


(103a 


(103t 
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In  line  182,  is  compared  to  10  .  If  d^  <  10  **r^  ,  then 

pq  -  pq 

little  confidence  can  be  placed  in  the  calculated  value  of  d2  because 
2  2 

r  and  t  are  very  close  to  each  other  in  (46b) .  If 
pq  o 


d 


2 


t 

o 


then  statement  52  stops  execution  because  an  accurate  value  of  1^  can 
not  be  obtained.  However,  if 


2  -5  2 

d^  <  10  r 

-  pq 


t  >±  A 
'  o'  —  2  q 


then  accurate  calculation  of  1^  may  be  possible.  If  |tQ|  -  y  is 
appreciably  greater  than  d^,  then  the  approximation  (90)  can  be  in¬ 
voked.  Line  184  puts  the  right-hand  side  of  (90)  in  W4.  Lines  187-190 
calculate  the  logarithm  term  in  (47a)  according  to  (91).  This  logarithm 
is  stored  in  W. 


Nested  DO  loops  45  and  46  calculate  (50)  for  n  =  Ml  +  M-l  where 
M  is  the  index  of  DO  loop  45.  The  index  K  of  DO  loop  46  obtains  &  in 
(50).  Nested  DO  loops  45  and  46  also  calculate  (91)-(93)  of  [1]  and  ( 91 )  —  ( 93 ) 
of  [1]  with  a  replaced  by  b.  The  results  of  these  calculations  are 
stored  according  to  (97)  and  (98).  With  regard  to  (68a),  lines  271-279 
put 


-  IT 


<t> 


(V 


1 


3  3 

2kpJ  J>= 

q 


V 


o  A  0 
q 


,3  3 
kpi 


log 


1  +  (V1-) 


2p  TT 

£_  \  L 


in  D8.  Line  280  puts  in  D9  the  term  which  is  present  in  in  case  4 
but  not  present  in  case  2.  Replacement  of  the  second  sum  with  respect 
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to  £'  in  (97)  of  [1]  by  the  exact  value  of  the  corresponding  integral 
with  respect  to  t'  and  use  of  (99)  of  [1]  give 


u<f>  (n.)  nt  (n  )  -jkR  0 

"  2  l  ko  ^  cos  4>a  costnf^)  l  A0,1  e  p 


5a 


£=1 


£’=1 


V 


kR 


+  D9 


p£'£ 


where 


(104) 


D9  =  - 


_2tt_  >  /V, 

*  08 


+  \  1  +  ('oil  )2 


2pq*l 


2pq^£ 


2irp 


2ttp  o 
1  +  ^ 


2rrp 

+  — r-^  log 


A  /  A  0 

wt  +  V1  +  (2ip") 

q  q 


(105) 


Line  280  puts  the  right-hand  side  of  (105)  in  D9.  DO  loop  67  adds  D8 
to  and  D9  to  G^.  DO  loop  67  also  sets  G^^,  G^j  and  G  ^  equal  to 

zero.  These  four  G's  are  set  equal  to  zero  when  p=q  because  the  exact 
values  of  the  coefficients  by  which  they  are  multiplied  in  (24)  are  zero 
then. 


The  index  L  of  DO  loop  13  obtains  £'  in  (35)  and  in  (62)— (63)  of 

[1].  With  regard  to  (36),  DO  loop  17  puts  G(u,<j)  )  in  GA(K) .  With  regard 

K. 

~^^n£'£ 

to  (64)~(66)  of  [1],  DO  loop  17  puts  (e  p  ^^^Rp£'£^  ®D(K).  If 

(59)  is  not  true,  line  318  sends  execution  to  statement  51.  Lines  319-340 
calculate  the  terms  in  (57)  which  do  not  involve  G(u,^)  and  the  terms  in 


(94)  of  [1]  which  do  not  involve  (e 
lates 


-jkR. 


p£'£ 


-)/(kRp£,£)’ 


n,  (V  ,2 

v*  A£  h 

L  3  3 

£=1  kaj 


in  D6, 


DO  loop  62  accumu- 


t  A .^)C 


£=1 


,3  3 
k  a£ 


+ 


2ka„ 


y'K 

«.3  5 

8k 


) 


in  D7, 


and 


71 


?  r 


ffit'i.  ..  r 


A 


in  D9. 


Line  337  puts  the  coefficient  of  n  in  (57c)  in  D8.  Line  338  puts  the 

portion  of  (57a)  which  does  not  involve  G(u,<j>0)  in  D6.  Line  339  puts 

^  2 

in  D7  the  portion  of  (57b)  which  involves  neither  G(u,<|)^)  nor  (n  +1). 
Line  340  puts  in  D9  the  terms  in  (94)  of  [1]  which  do  not  involve 

(e  p  )/(kR  £,g).  DO  loop  32  accumulates  G^,  G^,  and  G^  of  (36)  in 
H1A,  H2A,  and  H3A,  respectively.  DO  loop  32  also  accumulates  G^,  G,.,  and 
Gg  of  (64)-(66)  of  [1]  in  H4A,  H5A,  and  H6A,  respectively.  The  index  M 
of  outer  DO  loop  30  obtains  n  =  M  +  Ml-1.  If  (59)  is  true,  lines  368-370 
add  to  H1A,  H2A,  and  H3A  the  terms  in  (57)  which  are  not  present  in  (36) 
and  line  371  adds  to  H5A  the  terms  in  (94)  of  [1]  which  are  not  present 
in  (65)  of  [1].  Lines  372-380  do  the  sums  with  respect  to  £'  in  (35). 
Lines  381-386  do  the  sums  with  respect  to  £'  in  (62)-(63)  of  [1] . 

Lines  389-509  use  the  previously  calculated  G’s  of  (97)  and  (98) 

to  obtain  the  contributions  (23)  and  (24)  to  the  elements  of  the  moment 

matrices  Y  and  the  contributions  (48)-(51)  of  [1]  to  the  moment  matrices 
n 

Z^.  The  index  M  of  DO  loop  31  obtains  n  =*  M  +  Ml-1.  Lines  437-438  put 
G^a  and  G^  of  (54)-(55)  of  [1]  in  H4A  and  H4B,  respectively.  The  sub¬ 
scripts  Kl,  K2,...K8  defined  in  lines  452-459  have  the  same  meaning  as 
in  Table  2  on  page  50  of  [1],  The  variables  UA,  UG,  UB,  UH,  and  UF 
incremented  for  p=q  in  lines  425-429  are  intended  for  Y(K1),  Y(K2),  Y(K3), 

Y(K4),  and  Y(K9),  respectively.  Y(K9)  is  reserved  for  (Y^)  of  (24d). 

n  pq 


. .  -<J'  ■>'  vlv  /.vSTS.  a..  .v’iXfci 
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00  t  listing  OF  The  SUBROUTINE  YZ 

002  THE  SUBROUTINE  YZ  CALLS  THE  FUNCTION  BLOG 

303  SUBROUTINE  YZ I M l  .M2 . NP.N°H I , NT . I N, RH,  ZH.  X. A.  XT,  AT, Y,  Z ) 

004  COMPLEX  Y<  1600)  ,  7. 1  t  600)  .  U.  U  t  .  U2  *  U3  .U4.U5  .U6  ,  HI  A  ,  H2  A,  H3  A.  G  A(48) 

0  05  COMPLEX  GBI48). GCI48 ),  GD <48 ) .GE ( 48) . H l B . H2B . H3B. H 1C, H2C. H3C. H4A 

006  COMPLEX  N5A.H6A ,H4B  ,HSB .H6B , UA.UB , UC . UD , UE.UF , Gl A 1 1 0 ) . G2 A ( 10 1 

3  07  COMPLEX  G3  A <  1 0  )  .  Gl  B  <  10  ).  G2 E( 10 ) , G3BI 10 ) . G  1C (  10 ) . G 2C(  1 0) .G3C{  10) 

0  08  COMPLEX  G4AI 10) .GSAI 10) t  G6AI 10) «G4B( 10  I , G5B( 1 0 ) . G6BI tO).CMPLX 

009  COMPLEX  UG.UH 

010  DIMENSION  RHI43)  .ZHI  43)  .  X< 48 ) , A{  48) ♦ XT < 1 0 ) . AT ( 1 0 ) » RS  (42  )  .  ZS  < 42  ) 

01  l  D I  MENS  ION  D<42) . DR(42),  DZC42),  DM (42) ,C 1( 48) .C2( 48).C3(48) .C4C200) 

012  DIMENSION  C 5( 20 0 > »C6( 200 ) . R2 ( 10 > . Z2 <  l  0 ) ,R7 ( 1 0 ) ♦ Z7 ( 10 ) 

0  t  3  CT=2. 

014  CP  =  .  I 

015  DO  10  I =2 ♦ NP 

016  I  2= (— l 

017  RS(  (2 )=.S*(RH{ I )+RH(  12) ) 

018  ZS< I 2)=.5*(ZH(l )+ZH( (2) ) 

0  19  Dt =.5*( RH< I )-RH< (2) ) 

0 20  D2=. 5*( ZH( I )-Zri ( (2) ) 

021  D(  (2  )=SQRT <01*D  1  +  02*02) 

022  DR(C2)=Ot 

023  DZ  (  12  )-02 

024  DM( I 2)=D (I 2)/«S{ 12) 

0  25  10  CONT  INUE 

026  M3=M2-Mt+l 

027  M4=Ml-l 

328  PI2=l. 570796 

029  PP=  9  *869604 

0  30  DO  1  1  K=  t  •  NFHt 

031  PH=P(2*(X(K)<  „) 

932  Cl(K)=PH 

033  C2(K)=PH*PH 

334  SN=S IN( .5*PH) 

0  35  C3(K  )=4.  *SN*SN 

336  A l =P (2*  A  (  K  ) 

037  D4=,5*A  l*C3(K) 

038  D3-At*COS( PH) 

039  D6=A  t*6  IN*{PH) 

040  M5=K 

041  DO  29  M=t,M3 

342  PHM=(M4+M)*PH 

043  42=CO$(PHM) 

044  C4(MS)=D4*A2 

045  C5( M  5  )=D  5*A  2 

046  C6(MS)  =  C6*S  IN(PHM) 

047  M5=M5«NPHt 

048  29  CONTINUE 

049  H  CONTINUE 
050  PNl=.7853982*IN 

051  PN2=8.*PNl 

052  U=(0.»1.) 

053  Ul=.5*U 

0  54  U2  -2  »  *U 

055  MP=NP— l 

056  MT=MP-l 

057  N=MT+MP 

0  58  N2N=MT*N 

059  N2-N*N 

060  JN=-t-N 


*+.*.  » 


"  4  :  **«<****«»* 


DO  t 5  JO  »  l  »MP 
KQ=2 

IF(JQ.EQ.l)  KQ=l 
(F(JQ.EQ.MP)  KQs 3 
Rl=RS(  JQ ) 

21 =ZS( JQ) 

D  t=D( JO} 

C2=0R(JQ  ) 

D3=DZ< JQ) 

D4=D2/Rl 
D5=D  t/Rt 
SV=02/0l 
C  V=0  3/D  l 
Pl=PNl*Dl 
P2=PN2*OE 
P3=2  ,*Dl 
P4=2.*04 
P5=D4*D4 
P6=0 I*D l 
P7=Pi*0l 
T6=GT*Dl 
T  62=  T  6+0 l 
T  6'2=T62*  T62 
R6=CP*Rl 
R62=R6*R6 
PlA=0  . 

pi e=o. 
pic=o  . 

00  12  L=l.NT 
D6=XT(L ) 

R2 (L) =Rl +D2*D6 
Z2(L )  =7 1  +D3*D6 
07=P1 *AT<L)/R2(L> 

D8  =  t  .-06 
09=08*07 
PIA  =  P  LA  +08*09 
06=1 .+06 
Pt B=o 18+06*09 
Pt  C  =  PIC  +  C6*06*D7 
12  CONTINUE 

00  16  I P= l . MP 
R  3  =R  S I  IP  ) 

Z3=ZS(  IP) 

R4=R  l— R 3 
Z4=Zl-Z3 
U3=D2*Ul 
U4  =03  *U  I 
00  40  L=  l  .NT 
07  =R2  <L)-R3 
08=Z2 (L) -Z3 
R7 (L  )=R3*R2 (L) 

27(0=0  7*07  +  08*08 
40  CONTINUE 

PH=R4*SV+Z4*CV 
A 1  =ABS(PH> 

A 2  =A B  $ (  R  4* C  V  “Z4*  S V) 
D6=A2 

l  F  ( A  t.LE.Ol)  GO  TO  26 
06= A l -0 l 

D6=SQRT( 06*06+ A2*A2) 


26  IF( IP  .ME  .J0.AN0.<R6.GT .D6.0R .T6.LE.D6) )  GO  TO  41 
Z5=R4*R4+Z4*Z4 
RS=R3*Rt 
PHM=.  5*R3*SV 
00  33  K=t.NPHt 
A  1  =C  3(  K) 

RR=Z5  +  R5  *Al 
HI A=0. 

H2A=0  . 

H3A=0. 

H4A=0 . 

H5A=0. 

IFIRR.LT  .T62)  GO  TO  34 
OC  35  L= l « NT 
W=Z7(t_)4-R7(t_)*Al 
R=SQRT(W) 

SN=-S  IN  <R ) 

C$=CGS(R) 

06=AT(L )/R 

Hie=D6/W*CMPUX(CS~R*SN.SN+R4CS ) 

H  l  A=H  18+H  tA 
H2E=XT(L  >*H1 B 
H2A=H2B+H2A 
H3 A=XT(L )*h28+H3A 
H48=D6*CMPLX(CS. SN> 

H4A=H48+H4A 
H5A=XT( L  )*H4B+H5A 
35  CONTINUE 
GO  TO  36 
34  DO  37  L=t.NT 

W=Z7(L>+R7(L)*Al 

P=SQRT(W) 

1F(R.GT..5)  GO  TO  14 

CS*R* (W* {  .6944444E-2-W*. 1 736 11 l E-3 5- . 1 25 ) 

SN=W*( .33333336- l-W*. t 190476E-2)-.  3333333 
HI 0=ATiL  )*CMPt_X  (CS.SN) 

CS=R*(W*( .4166667E- 1-. 1 3 88889E-2*W >-.5 > 

SN=W4(W* ( . 19841 26E-3*W-.8333333E-2)+. 1666667)- 1 . 
H4B=AT<  L)*CMPt_X(CS»SM 
GO  TO  43 
14  SN=-SIN(R) 

CS=COS(R ) 

06  =AT (Li/R 

H 18=06* ( (CMPLX(CS-R*SN,SN+R*CS)-l .)/W-.3) 
H4B=D6*CMPLX(CS-  l  ..SN) 

43  HI  A=H  184-Hl  A 
H2B=XT(L )  *H  l  B 
H2A=H2B+H2 A 
H3A=XT(L )*H28+H3A 
H4A=H4B+H4A 
H5A-1  XT  (  L  )*H48+N5A 
37  CONTI NUE 

A  l  =PH4PHM  *A 1 
A2=AFJS(  At ) 

R=RR-A2* A2 
06=A2— 01 
D7=A240 1 
062=06*06 
072=07*07 
08  =  SQRT ( 062  4R) 


on;-*- 

..  ..  .  ..  •*  1  i.'  .-.i-  . 


v 
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181 
l  82 
183 
l  84 
185 
t  86 
187 
1  88 
1  89 

190 

191 

192 

193 

194 

195 

196 

197 

198 

1  99 
200 
201 
202 

203 

204 

205 

2  06 
2  07 
208 
209 
2  10 
2  l  1 
212 
2  13 
2  14 
215 
2  16 
2  17 
2  18 
219 
2  20 
221 
2  22 

223 

224 
2  25 
226 

227 

228 
2  29 

230 

231 

232 
2  33 
2  34 
2  3S 
2  36 
2  37 
2  38 
2  39 
2  40 


D9=SCRT (072+R) 
lF(R-lRR*t.F-5))  52.52,53 

52  IF1D6.LT.0.)  STOP 
Vl4  =  .5/D62-.5/D72 
GO  TO  54 

53  W4=( 07/09-06/08 >/R 

54  tF(D6.GE.O.)  GO  TO  38 
Vi=ALGG(  (  07+09)  *< -D6  +  08  )/&> 

GO  TO  39 

38  W=ALC!G(  ( 07  +  09 )/{  06+08 )  ) 

39  Wt=(W4-fr.5*W)/Ol 
W5= A2/01 

W 2=1  .5*109-08)- l  ./09+1 ./08)/P6-W5*Wl 

W3  =  (  .25*107*09-06*08  )  +W-R*  <W 4+  .25*W  )  )/ P7-W 5* ( 2 .  *W24-W5*W  l  ) 
W4=W/Dl 

W5=l 09-08- A2*W)/P6 
IF1AI.GE.0.)  GO  TO  27 
W2=-tf2 
V»5*-  W5 

27  H1A=VH+HIA 
H2A=VK2*H2A 
H3A=W3+H3A 
H4A  =  \»4*H4A 
H5 A=W5+H5A 
36  GA 1 K ) =K l A 
GB(K)=H2A 
6C(  K )=H3A 
GD(K)=H4 A 
GE(K) =H5 A 
33  CONT  INUE 
Kt=0 

00  45  M=l,M3 
HI A=0. 

H2A=0  . 

H3A=0. 

HI  0=0  . 

H2B=0. 

H3  0=O  . 

HlC=0. 

H2C=0  . 

H3C=0. 

F4 A=0  . 

H5A=0. 

H6A=0  . 

H48=0. 

H58=0  . 

H6B=0. 

DO  46  K=  l.NPHt 
Kt  =Kt  +  l 
06=C4<K  l  ) 

07=C51IU  ) 

08=C6tKt  ) 

U A=G  # ( K ) 

UB=GB(K) 

UC  =  GCOO 
UO=GO 1 K ) 

Ue=GSlK> 

H l A  =  0  6*U  A+H  IA 
H2 A=C7*U A  +  H2A 
H3A=08*UA*H3A 
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241  HIB=D6*UB4-HIB 

242  H2B=074Ue+H2B 

243  H3B=Da*US+H38 

244  HIC=06*UC+H  1C 

245  H2C=D7*UC+H2C 

246  H3C=08*UC+H3C 

247  H4 A=D6*U0+H4A 

248  HSA=07*UD+HSA 

249  H6  A=D8*  U0+H6A 

250  H4B=D6*UE+H4B 

251  H58=D7*UE+H5B 

252  H6B=D6*UE+H6B 

253  46  CONTINUE 

254  GIA(M)  =  H1A 

255  G2A(M)=K2A 

256  G3A(M)=H3A 

257  G1 8 ( M ) =H1 B 

258  G2B( M )=H2S 

259  G3B(M>=H3B 

260  G IC(  M  )  =  H  1C 

261  G2C(M)=H2C 

262  G3C(M)=H3C 

263  G4A(M)=H4A 

264  G5A(M)=HSA 

265  G6A(M)=H6A 

266  G4B(M)=H4B 

267  G58(M)=HSB 

268  G6B  C  M)  =H6B 

269  45  CONTINUE 

270  IF(IP.NE.JQ)  GO  TO  47 

271  A I =05  ADS 

272  D8=0. 

273  D9=0  . 

274  DC  63  K=l,NPHt 

275  D8=D8+AIK ) /SORT ( C2< K) 4A U 

276  D9=D9+ A  (  K )  *  ei_OG  (05/ Cl  (  K|  ) 

277  63  CONTINUE 

278  A2=3, 141 593/05 

279  D8=(BLOG(A2)-PI2*C8>/(R5*R  I) 

2  80  09=2. /Rl*( BLOG ( A2)  A2*  BLOG  ( l  ./A2 ) )-3 .1 4 1593/DI*D9 

28  1  DO  67  M=  I.M3 

282  Gl A( W)=De+Gl ACM) 

283  G 2  A  I  M  )=0  • 

284  G2e(M)=0. 

285  G2C(M)=0* 

286  G3A(M)=0. 

287  G5A(M)=DS4G5A(M) 

288  67  CONTINUE 

289  GO  TO  47 

290  41  DO  25  M=t ,M3 

291  GtA(M)=0. 

292  G2  A ( M ) =0  « 

293  G3AI M  )=0  • 

294  GIB(M)=Q. 

295  G 2BI M )  =  0  » 

296  G3 E I M ) =0  « 

297  G1C(M)=0. 

2  98  G2  C  ( M )  =0  » 

299  G3C(M>=0. 

300  G4A(K)=0. 


G5A(  M )=0  . 

G6  A ( M ) =0  . 

G4B( M  )=0  • 

G5e(M)=0  . 

G68(  M )=0  . 

25  CONTINUE 

DO  1  3  L=  l,NT 
R5=P7(L) 

ZS=Z7(L> 

DO  17  K= l <  NPHT 
W=ZS+R5*C3<  K) 

R=$QRT(W  > 

SN=-SIN<  R) 

CS=COS(R  ) 

GAIK )=CMPLX(CS-R*SN,SN+R*CS> /{W*R} 

GD(K)=CMPLX (CS.SN J/R 
17  CONTINUE 

lF(R62.t_E.Z5)  GO  TO  51 
06=0. 

07=0  . 

09=0. 

DO  62  K=  l.NPHt 
W  2  =C  2  (  X  ) 

W= 1 . /<Z5+R5*W2> 

Wl=A(K)*SQRT(W> 

0  6=0 6 +W  l  *W2*W 

07=0  7FW  I  *(  .  S+W*  (l.-*-.  125*  W*R5*W2*W2>) 

09=09 +W  l 
62  CONTINUE 
W l=R5/Z5 
„2=PP*VJl 
W=  SORT ( W  2 ) 

Vi3=t  .  +W2 
R=SQRTIW3) 

Vi  4  =S  CRT  (  R5  I 
W  5= A LOG ( VtR  ) 

08=-P 12*  C6- ( W/R-W5 )/ ( R5*W4  ) 

06= .5*0 e 

07  =  1 <W/R*IWl-( .  125*.  1 666667**2  >/W3>-*-. t  25*W  5)/R  5+ .  5*W  5  > /W4-P I  2*D7 
09  =  WS/W4-Pt  2*09 
51  Al  =  AT  (L  ) 

A 2  =XT ( L  MAI 
A3=XTCU)*A2 
K  1=0 

DO  30  M=1 ,M3 
W=M+M4 
HI A=0  . 

H2A=0  . 

H3  A=0  . 

H4A  =  0  . 

H5  A=0  . 

H  6  A=  0  • 

DO  32  K=l,NPHt 
Ki  =  KlH 
HI B=GA( K ) 

W4=C  4 ( K l ) 

N5=C5 (Kl > 

V»6=C  6 ( K I  ) 

Ht  A=*4*H  10  +  HlA 
H2A=W5*H IB+H2A 


H3  A»W6*W IB  +  H3A 
Hi B=GO< K > 

H4  A=W4*Hl  B*H4A 
H5A-W5*H IB+H5A 
H6  A  =  *J6  *H  1 8  +  H6A 
32  CONTINUE 

IF ( R62  c LE • 15 }  CO  TO  44 
Ht A=D6+H  t A 

H2A=07-(W*W+t  .)  *i5S*H2A 

H3A=W*D8*H3A 

H5 A=D9 +H5A 

44  G1A(M|=AI*HIafGIA(M) 

G2A(M)=A1*H2A«-G2A(M> 

G3A( M)=A l*H3A+G3A{Ml 
G18<M)=A2*H1  A+G1BIM) 

G2B ( M)=A2TH2AvG26  <MJ 
G3B  (  M  >-A2*H3A«-G3B<M) 

GIC<  M)=A3*H l A*Gl C(M» 

G2C( M)= A?*H2A+G2C(M> 

G3CC  K)=A3*H3A*G3C(M> 

G4A (M )= A  l*K4A+G4A(M) 

G5A<  M)  =  Al*H5A*G5A(M» 

G6A  (  M  )=Al*H€A+G6A(M> 

G481 V)=A2*H4A4G4B(M) 

GSB  (M  )=  A2*H£A-i-G5B<M  ) 

G6BC  M  )  =  A2*H6  A4-G6  8  t  M  ) 

30  CONT  1NUE 
13  CONTINUE 
47  A2=0{  IP) 

A3=.5*A2 

W l=A  3*(  R  4*0  3— Z 4 *02) 

V»2=-A3*P3*D3 
A3=DZ(  IP) 

06=DF< IP) 

D7=Z  4*06 
09  =  D3*D6 

H  t  C=( R 1 *D9— 0  2* (R3*A3+07) )*U 

08  =A2*Dl 

H3C=08*U 

H2  C=Z4*  H3C 

H3C=D3*H3C 

W3=°3*( R4*A3— 07 ) 

W4  =P3  * (  C2* A3— 09  ) 

W5=P3*R  l* A 3 
AC  =0R(  IP  ) 

U5=A l*U3 
U6  =  .”*U4 
06=-02*A2 
0  7=D  I.  *  A  l 
A  3  =0  M ( IP) 

JM=JN 

00  31  M=UM3 
H2  As  G2  A  (  M  ) 

HIA=GIA(M> 

H2B=G2B ( M ) 

H  l  B=6  l  3  (  W) 

UC=W  t  *H2 A  +  W2*Hl A 
JB=W l*H28*W?*Hl 8 

UF=W3*(  H2A  +  D4*H2B)*W4*(  H?BT04*G2C( M) )*W5*( H 1  A* P4*H l B+P5* G ICC M) J 
UA=UC-U8 
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42t  UB=UC+U0 

422  UG=UA 

423  UH=UB 

424  IF(IP.NE.JC)  GO  TO  48 

425  UA=PtA<-UA 

426  UG  =  P  t  BFUG 

427  UB=PtS+Ue 

428  UH  =P  IC-HJH 

429  UF=P?+UF 

430  48  H3A=G3A(M) 

431  H36=  G3B ( M ) 

432  UC=H1C*H3A 

433  UD=HIC*H3B 

4  34  UE=H2C*  <H3A«-04*H3B) +  H3C* ( H3B +04  * G3C (M ) ) 


435 

H5A=G5A(M) 

4  36 

H5B=G58<  M) 

437 

H4A=G4A(MM-H5A 

4  38 

H4B  =  G4B<  M) +H58 

4  39 

H6A=  G6A (M ) 

440 

H6E=G68(  M) 

441 

H3A=U5*H5AFU£*H4A 

442 

HI B=US*H5B*U6*H4B 

443 

Ht  A=H3A~HIB 

444 

H2A=H3A+Hl B 

445 

H3A=-UI.*H4A 

446 

HI B=D6*H6A 

447 

W=H+M4 

481 

448 

A  l  A3 

482 

4  40 

H2B=D6*H6B— A l*H4A 

4  83 

450 

H3B=C7* ( H6A+044H68) 

484 

4  5  t 

H4A=W*D5*H4A 

4  85 

452 

KS  =IP  +  JM 

4  86 

453 

K2=Kl+t 

4  07 

4  54 

K3=Kl+N 

4  88 

4  55 

K4=K2+N 

4  89 

4  56 

K5=K2+MT 

4  90 

457 

K6=K4*MT 

491 

458 

K7=K3+N2N 

492 

459 

K8-K4FN2N 

4  93 

4  60 

<9=K8+MT 

494 

461 

GO  TO  {  18,20,19)  <  KQ 

495 

462 

18 

Y  (  K6  )=UC+UO 

4  96 

4  63 

ZI K6) =H  t  BFH2B 

497 

464 

IF t  IP.EQ.l )  GO  TC  21 

4  98 

4  65 

Y(  K3)=Y(K3)+UB 

499 

466 

Y  ( i<7  l=Y  (  K7  )  +UE 

500 

4  67 

Z( K3) =Z ( K3) FH2A-H3A 

301 

4  68 

Z( K7  )=Z(K7 ) +H3B— H4A 

502 

469 

IP.EO.MP)  GO  TC  22 

503 

4  70 

2  l 

Y I K  4 )=UH 

5  04 

4  71 

V { K8)=UE 

505 

472 

Z(K4)=H2  AFH3  A 

506 

4  73 

Z( K8)=H39+H4A 

507 

4  74 

GO  TO  22 

503 

4  75 

1 9 

Y(  X5) =Y( K5)+UC“UD 

509 

4  76 

Z( K5  )  =  Z(K5) FHIB-H2B 

5  10 

4  77 

IF(  IP.EQ.l  )  GO  TO  23 

5  l  l 

478 

Y(Kl)=Y(Kt H-UA 

5  t  2 

4  79 

Y(K7)=Y(K7 )+UE 

5  13 

4  BO 

Z ( K l )=  Z IK l ) FH  t  AFH3A 

514 

Zt  K7) =Z( K71+H3B-H4A 
IF (  IP  ,EQ  .MP  )  GO  TO  22 

23  Y<  K2) =  Y( K2J+UG 
Y(K8>=UE 

Z<  K2  )  =Z<  K2  ) +HI A— H3A 
Z(K8 )=H3B+H4A 
GO  TO  22 

20  Y(K5)=Y(K5)FUC-UD 
Y ( K6)=UC*U0 
Z(K5)=Z(K5) +HIB-H2B 
Zt  K6)  =Hl  84H2 B 
IF(  IP.EQ.l)  GO  TO  24 
Y  (  Kl  )  — Y  (  K l  >  -HJA 
Y (  K3)=M  K3)*UB 
Y(K7)=Y(K7)+UE 
Z  (  K  1  )=?.  \  K  l  )  FHIA+H3A 
Z(K3>  =  Z(K3)4F2  A— H3 A 
Z( K7)=Z< K7J+H3B-H4A 
IF(IP.EO.MP)  GO  TO  22 

24  Y( K2 )=Y( K2 ) +UG 
Y I K4 } =UH 
Y(Ke)=UE 

Z(K2)=Z(K2)+Ht A-H3A 
Z( K4  )=H2A+H3A 
Z(K8) =H3  B  +  H4 A 
22  V(K9)=UF 

Z< K9) =U2*(D8*( K5A404  4H5B  )- At  *H4A  ? 
JVi  =  JM+N2 
31  CONTINUE 
16  CONTINUE 
JN=JN+N 
t5  CONTINUE 
RETURN 
END 
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VIII.  THE  MAIN  PROGRAM  FOR  BOTH  THE  H-FIELD  ARP  E-FIELD  SOLUTIONS 

The  main  program  for  both  the  H-field  and  E-field  solutions 

obtains  the  present  H-field  solution  and  the  E-field  solution  of  [1] 
for  the  electric  current  on  a  conducting  body  of  revolution  immersed 
in  an  incident  plane  wave.  Input  data  are  read  from  punched  cards. 
These  input  data  are  the  same  as  those  for  the  main  program  for  the 
H-field  solution  which  is  described  and  listed  in  Section  VI  of 
Part  Two. 

The  main  program  for  both  the  H-field  and  E-field  solutions 

calls  the  subroutines  YZ,  PLANE,  DECOMP,  and  SOLVE.  The  function 

subprogram  BLOG  is  also  needed  because  it  is  called  by  the  sub¬ 
routine  YZ 


where 


Minimum  allocations  are  given  by 

COMPLEX  Y(N*N) ,  Z(N*N),  RE(2*N) ,  R(2*N),  B(N) ,  C(N) 
DIMENSION  XT (NT) ,  AT (NT) ,  X(NPHI) ,  A(NPHI) ,  RH(NP) , 
ZH(NP) ,  IPS(N) ,  IPT(N) 

N  =  2*NP-3 


The  t  and  cj)  components  of  the  present  H-field  solution  for  the 
electric  current  are  calculated  from  (84)  and  (85)  .  The  t  and  <t>  com¬ 
ponents  of  the  E-field  solution  of  [1]  for  the  electric  current  are 
also  calculated  from  (84)  and  (85).  In  the  present  H-field  solution, 
the  coefficients  1^  and  1^  in  these  equations  are  the  pth  elements  of 
the  vectors  1^  and  1^  which  satisfy  the  n=l  equation  in  (10).  In 
the  E-field  solution  of  [1],  these  coefficients  are  the  elements  of  the 

vectors  1^  and  1^  which  satisfy  the  n=l  equation  in  (6)  of  [1] .  DO  loop 

28  prepares  RH  and  ZH  for  use  in  the  subroutines  YZ  and  PLANE  by  multi¬ 
plying  them  by  k.  With  regard  to  (10),  line  42  puts  Y^  of  (86)  in  Y. 

With  regard  to  (6)  of  [1],  line  42  puts  Z^  of  (96)  in  Z.  Line  45  calcu¬ 
lates  IPS  and  changes  Y.  Line  46  calculates  IPT  and  changes  Z.  With 
regard  to  (10),  line  47  puts  the  excitation  vectors  1^  and  1^  for  the  0 
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polarized  incident  plane  wave  (69)  in  RE.  With  regard  to  (6/  of  [1] , 
line  47  puts  the  vectors  and  -V^  for  the  0  polarized  incident  plane 
wave  iaR.  Line  47  also  stores  vectors  for  the  cj>  polarized  incident  plane 
wave  (70)  further  on  in  RE  and  R  but  these  vectors  are  not  used  in  the 
main  program. 

In  DO  loop  36,  JHE  =  1  obtains  the  present  H-field  solution  for 
the  electric  current  and  JHE  =  2  obtains  the  E-field  solution  of  [1] 
for  the  electric  current.  In  line  55,  the  output  IPS  and  Y  from  the 
subroutine  DECOMP  is  fed  along  with  N  and  RE  into  the  subroutine  SOLVE. 
SOLVE  puts  the  solution  1^  and  1^  to  (10)  in  C.  Lines  59-64  put  V^ 
and  V^  in  B.  In  line  67,  the  output  IPT  and  Z  from  the  subroutine  DECOMP 
is  fed  along  with  N  and  B  into  the  subroutine  SOLVE.  SOLVE  puts  the 
solution  1^  and  1^  to  (6)  of  [1]  in  C.  The  t  and  <j>  components  of  the 
normalized  electric  current  are  printed  out  under  the  headings  JT  and  JP, 
respectively.  The  normalization  is  the  same  as  in  the  main  program  for 
the  H-field  solution  which  is  described  in  Section  VI  of  Part  Two.  The 
sample  input  and  output  data  are  for  the  sphere  examples  of  Figs.  2 
and  4. 
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OOIC  LISTING  OF  THE  MAIN  PROGRAM  FOR  BOTH  THE  H-FtELO  AND  E-FtELD  SOLUT  tONS 

002  C  THE  SUBPROGRAMS  YZ. BLOG, PL ANE, DE COMP ,  AND  SOLVE  ARE  NEEDED 

00 3/ /PGM  JOB  ( XXXX.XXXX, l , t)  ,‘MAUTZ  .JOE* ,REGt 0N=2OOK 

004//  EXEC  WATFtV 

005 /  /GO  •  SYS  IN  DO  * 

006S JOB  MAUTZ.T  IME= 5, PAGES3 60 

007  COMPLEX  Y{ 1600) ,Z{ t600) .RE (240) ,R(240) .8(40) ,C(40 J.U.Cl 

0  08  Dt ME NS  ION  XT ( 10 ) , AT ( 10 ). X ( 48 ) . A( 48),RH( 43) .ZH( 43) , THR( 3)  . IPS( 40) 

009  D ( MENS ( ON  tPT(40) 

0 10  READ (1.15)  NT.NPHt 

011  15  FORM  AT(  2  t  3) 

012  WRITE{3,  30  )  NT.NPHt 

0t3  30  FORMAT ( •  NT  NPH t • /t X , t3 , 15 ) 

014  PE AD (1*10) (XT ( K ) ,K=  t . NT) 

0  t  5  READ( t.tO) (AT(K)  .K=l  .NT) 

0  16  10  FORMAT (5E14 .7) 

017  WRITE <3. tt) (XT(K) »  K  =  l . NT ) 

018  WP.  (TE(3,  12)(AT(K)  »K-  1  .NT  ) 

019  ll  F  O  RMAT  (  *  XT  •  /  (  t  X  »  5  E  t  4  »'  7  )  ) 

020  12  FORMAT ( *  AT • /( l X , 5E l 4 . 7) ) 

021  REAOC  l.  10) (X (K)  »  K=1  «  NPH  1 ) 

022  READ (tjtO)(A(K)*K=l, NPH l ) 

023  WR  (  TE  (  3  »  1 3  ) (X(K) . K= 1 . NPH l ) 

024  WRITE  (3,  1 4  )  {  At  K  )  ,  K=  1 ,  NPH  I  ) 

025  13  F ORMAT ( *  X* /(IX .5E14.7)) 

026  14  FORMAT ( •  A« / ( 1 X. 5Et 4 . 7 ) ) 

027  READ (1.16)  NP,BK,THR(l> 

028  16  FORMAT (13.2E14.7) 

029  WR  ITE(3 , 17 )  NP.BK.THR(l) 

030  17  FOPMAT ( *  NP« , 6 X, • BK * , 12 X . *THR • / l X . 1 3. 2E l 4 . 7) 

031  REAO(l.l8)(RH(t).t=l.NP) 

032  READ( l,  t8)(ZH( t ).t=l.NP) 

) 33  18  FO  RM  AT ( 1 0F8 .4 ) 

034  WR  [TE(3.  19)(RH(  f  )  ,  [  —  l  , NP ) 

035  WR t TE (3 ,20) ( ZH(  l),t  =  l»NP) 

036  19  FORMATt*  RH • /(  IX. 10F  8,4)) 

037  20  FOPMAT ('  ZH • / ( t X , 10 F8  .4 ) ) 

038  DO  26  J=l.NP 

039  RH  <J)=BK4RH( J) 

040  ZH( J)=8K*ZH( J> 

041  28  CONTINUE 

042  CALL  YZ(  1.1  .NP.NPHl  .NT , l ,RH. ZH , X , A , XT, AT , Y,  Z  ) 

0  43  MT =NP— 2 

044  N=2*MT+l 

045  CALL  DECOMP (N, IPS.Y) 

046  CALL  DEC  OMP(N.IPT.Z) 

047  CALL  PLANE(  l.l,  l.NP,  3.  NT.PH,  ZH,  XT,  AT.THR.RE.R) 

048  U=(0..l.) 

049  DO  36  JHE= 1.2 

050  GO  TO  (33,34) » JHE 

051  33  WR  (TE(3.29)  (  Y(  J  )  ,  J=  UN) 

052  29  FORMA  T (  '  Y » / ( 1 X , 6E 1 l . 4 ) ) 

053  WR  (TE  (3,35)  (PE(  J  ).J=  l.N) 

054  35  F ORMA T ( *  R« / ( l X , 6EI t  .4  )  ) 

0  55  CALL  SOLVE  (N.tPS.Y.RE.C) 

056  GO  TO  32 

057  34  WP(TE(3,3l  >(Z(J  >  »J=1»N) 

058  31  F  ORMAT ( 1  Z • / ( 1 X . 6EI 1 .4 ) ) 

059  DO  22  J=l,MT 

060  B( J) sR( J) 
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061 

J  1  =  J  +  MT 

062 

B(  Jl)=-R(Jt) 

0  63 

22 

CONT  tNUE 

06* 

BIN)  =-R(N) 

0  65 

WRITE  (3. 23  )<B<  J  ),  J=  t  ,  N  > 

066 

23 

F  ORMAT (  '  F  * / ( IX . 6E1 1  ,4)) 

0  67 

CALL  SOL V£(N, IPT.Z.B.C) 

0  68 

32 

WR I TE (3,21) 

0  69 

21 

FORM  AT ( •  REAL  JT  {MAG 

JT 

070 

DO  24  J=  1  »  MT 

071 

C 1=2  • /RH ( J  +  1)*C(J) 

0  72 

C2=CABS<Cl ) 

0  73 

WR  IT  E  (  3.  25  )  C1.C2 

0  74 

25 

F  CKMAT ( 1 X.3E11 . 4) 

075 

24 

CONT INUE 

9  76 

WR(TE(3 .26) 

077 

26 

FORMAT ( •  REAL  JP  (MAG 

JP 

0  78 

MP=NP-l 

079 

DO'  27  J=  l.MP 

080 

Ct=4./(RK( J)  +  RH (J+l ) >*U*  C< J  +  MT ) 

081 

C2=CAeS(Cl) 

082 

WR IT  E (3 , 25 )  C1.C2 

083 

27 

CONT  CNUE 

0  84 

36 

CONT  tNUE 

0  85 

STOP 

0  86 

ENO 

S  DATA 

MAG  J  T*  ) 


MAG  JP') 


2  20 

-  0  .5773503E+00  0  .5773503E+00 
0.  1000000E+01  0  <  1  000000E+0 1 

-  0.0031 2 86 E +00—  3  .06307 1 9E+O0-O . 0 1 22344E+00-0 . 8391 t70E+00-0  .746 3319  E+00 

-  0. 636  05 37 E *00-0 .5108670E+00-0  .37370 6 1 E+OO-O . 2277S59E+00-0 . 7652652E-01 
0.76526E2E-01  0 . 2277859E  +  00  0 • 373706 1 E+00  0 . 5  1 0867OE+0O  0  .6360537E+00 

0  .9  122344E+00 
0  o  6  267205E— 0 1 
0 . 1  42096 l E+00 
0. I  423961 E+00 
0  .6  267205E—  0 1 


0 . 746331 9  E  +00 
0. 1761401E-0  t 
0.1  l 9 l 945  E +00 
0. 1527S34E  +  00 
0 . 1019301  E+0  0 


0 .8391 l 70E+00 
0.40601 43E—  0 1 
0 . 131 6  8  86E+00 
0 , 149173 OE +00 
0 .8327674E-01 


0 .96397 19E+00 
0. 8327674E-0I 
0 .  1491730E+00 
3. 13168 86 E+00 
0.4060 143E-01 


0 . 993 1 286E+  00 
0  .1019301 E+00 
0. 1 5?  7534E+00 
0 . 1 181945E+O0 
0. 176 140  IE- 01 


16  0.  1256637E+01  0.3141 593E+0 1 


0.9511  0.9945  0.9945 


0.3090  -0. 1045  0. 1045 


0  .0000  0  .2079  0.4067  0.5878  0.  743t  0.8660 

0.8660  0.7431  0.5878  0.4067  0.2079  0.0000 

-1  .0000  -0.9781  -0.9135  -0.8090  -0.669t  -0.5000 
0.5000  0.6691  0.8090  0.91  35  0.9871  V  .0000 

S  STOP 
/  * 

// 

PRINTED  OUTPUT  CONSISTS  OF  THAT  FROM  THE  MAIN  PROGRAM  FOR  THE 

SOLUTION  IN  SECTION  Vt  OF  PART  TWO  PLUS  THE  FOLLOWING  LINES. 


H-F  (ELD 
Z 

0 .6733E-0  1-0  .  U0  8E+0  2 
0.4415E-01  0.311  2E+0  0 
0 . 1437E-0 l 

-  0 . 27  39E-0  2 

-  0 .7250E-0  2 
0. 27( 5E+01 


0.6205E-0  l 
0. 3365E-0 l 
0  •  6434E-0  1  0.6919E-02 
0 . 3823E— 0 1— 0.5271E-02 
0 .3855 E- 01— 0.758 4E -02 
0 . 6766E— 0  l  0.4763E+00 


0.31  62E+01 
0 . 1 5  93E+O0 
0-4F79E-01 
0.3745  E—  0 1  - 


0 . 5399E- 0  l 
0  .2345E-0  l 
0. 1242E-02 
0  .6728E-02 


0.3880E-01-0. 1503E+02 


0  .7427E-0  l 
0. 4254E-0 l 
0 .40  74E-0 l 
0,  3996E-0 1 


0.531  5E—  0  l 
0. 331 7E-01 
0 . 17O7E-0  1 
0 .8753E-02 


C.5232E-0  l 
0.4167E-0  l 
0.4026E-0  1 
0.3906E-0  1 


0  .64 12E-01 
0.46  56E-01 
0 .2706E-01 
0. 1339E-01 
0.7696E-02 


0 .1469E+00 
0.450 l E— 0 1 
0  .4120E-0  1 
0 » 3982E— 0 1 


0. 8394E+00 
0.9534E-01 
0. 4132E-QI 
0.3773E-01 
0.6952E-0I 
0.5914E-01 
0 • 3977E-01 
0.2165E-01 
0. 1064E-01 


0.951  l 
0.3090 
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- 

0 . 2754E+00-0 

.  7539E+00— 0 

.  3 1 20E+00— 0 .68  It  E+OO-O  .3507E+00-0 . 5623E+00 

- 

0. 3645 E +00-0 

«  4085E+00-0 

.3286E+00-0. 2433E+00- 0.231 3E+00-0. l0t2E+00 

- 

0.83SIE-01-0 

.  1840E-0 t 

0 

•  83 51E-01  — 0.1840  E-C  t  0  .23 l 3E+00-0 . tO 1 2E+00 

0  .3286E+00-0 

. 2433E+0  0 

0 

.3645E+00-0.4085E+00  0 . 3507E+ 00-0. 5623E+00 

0.3106E+00-9 

.681 7E+00 

0 

.  2697E +00-0. 75S9E+00  0. 780 l E+00-0. 2694E+00 

0.7641E+09-0 

.  3t  l  tE+00 

0 

.  7268E+00— 0.3900 E+00  0  .6583E+00-0 . 4958E+00 

0  .550  5E -^10-0 

.612  8E+00 

0 

.  3995E  +00-0.7201E  +  00  0 . 21 07E+00-0. 7956E+00 

- 

O.OOOOE+OO-O 

.  8227E+00— 0 

•  210 7E +00— 0.79 56 E+00-0  .39 95 E+0 0-0 .720  t E+00 

- 

0.5S05E+00-0 

*  6  1 28E+0  0- 0 

•  6583E +00-0. 49  58E+  00-0. 7268E+00-0. 3900E+00 

- 

0  •  7764E+00— 0 

.  3 1 1 6E+00-0 

.  7 7 93 E+0 0—0 .26  26  E+00 

REAL  JT 

tM AG  JT 

MAG  JT 

- 

0. 8238E+00-0 

.  1 97 1 E+0  t 

0 

.  2 1 37E+0 1 

- 

0.9777E+00-0 

. 1907E+O t 

0 

.2143E+01 

- 

0. 1 1 86E+01-0 

. 1802E+0 1 

0 

•  2  1 58E+0 1 

- 

0. 1421E+0  t-0 

. 1620E+0 1 

0 

. 2  1 54E  +0  l 

- 

0.  1629E+01-0 

. 1337E+0  l 

0 

.2 l 08E+0 l 

- 

0 . 175  IE+0  t-0 

. 9546E+00 

0 

. 1994E+0 l 

- 

0. 1733E+01-0 

. 50 l 3E+00 

0 

.1804E+0  t 

- 

0. 1545E+0  t-0 

. 2843E-01 

0 

.  1545E  +  01 

- 

0.  U96E+01  0 

. 404 5 E+00 

0 

•1263E+0 l 

- 

0.7315E+00  0 

. 7  493E+00 

0 

.  1047E  +  0  l 

- 

0*21 98E+00  0 

•  9839E+00 

0 

.  1008E  +0  t 

0  »  26  53E+00  0 

. tt  t  6E+0  l 

0 

. t 147E  +  0  l 

0.66  3 1 E+00  0 

. 1182E+01 

0 

,1 355E+0  1 

0.9285E+00  0 

. 1 196E+0  l 

0 

.  t5l4E+0t 

PEAL  JP 

IMAG  JP 

MAG  JP 

0 ♦ 7997E+00  0 

. 1970E+0  l 

0 

.  2  1 26E  +  0 t 

0.9274E+00  0 

-  1833E+0  l 

0 

.2054E+01 

0  •  tO 7 IE+0 l  0 

. 1645E+01 

0 

.  1 963E+0 l 

0  .122 IE+0  l  0 

. 1362E+0 1 

0 

.1829E+01 

0 . 13 1  2E  +0  t  0 

. tOt  4E+0  t 

0 

.  l658E  +  0t 

0.1297E+01  0 

.6571 E+00 

0 

.  1454E+0  t 

0 . 1 1 66E  +  0 1  0 

, 3656E+00 

0 

. 1222E+01 

0  .956  l  E+00  0 

.20585+00 

0 

.9780E+00 

0.  7409E  +  00  0 

,  2063E+00 

0 

.  769tE+00 

0. 592 t E+00  0 

• 3454E+00 

0 

.6855E +00 

O.SEOeE+OO  0 

.  5635E+00 

0 

•  7880E  +00 

0.61  20E+00  0 

•7924E+00 

0 

.  100  IE  +0  t 

0. 7237E  +  00  0 

.  9694E+00 

0 

«  l  2  t  OE+Ot 

0.8687E+00  0 

. It  13E+0  l 

0 

. 1412E+01 

0.9847E+00  0 

.  1 1 94E+0 t 

0 

. l547E+0t 
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